Numerical studies of the phase diagram of layered type II superconductors in a 
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We report on simulations of layered superconductors using the Lawrence-Doniach model in the 
framework of the lowest Landau level approximation. We find a first order phase transition with 
a B(T) dependence which agrees very well with the experimental "melting" line in YBa2Cu307_a. 
The transition is not associated with vortex lattice melting, but separates two vortex liquid states 
characterised by different degrees of short-range crystalline order and different length scales of cor- 
relations between vortices in different layers. The transition line ends at a critical end-point at 
low fields. We find the magnetization discontinuity and the location of the lower critical magnetic 
field to be in good agreement with experiments in YBa2Cu307_a. Length scales of order parameter 
correlations parallel and perpendicular to the magnetic field increase exponentially as 1/T at low 
temperatures. The dominant relaxation time scales grow roughly exponentially with these correla- 
tion lengths. We find that the first order phase transition persists in the presence of weak random 
point disorder but can be suppressed entirely by strong disorder. No vortex glass or Bragg glass 
state is found in the presence of disorder. The consistency of our numerical results with various 
experimental features in YBa2Cu307_^, including the dependence on anisotropy, and the temper- 
ature dependence of the structure factor at the Bragg peaks in neutron scattering experiments is 
demonstrated. 



PACS numbers: 74.20.De, 74.25.Dw, 74.25.Ha 
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I. INTRODUCTION 



In the mean-field limit the phase diagram of type II 
superconductors has two phases: the normal state and 
the mixed state in which the lines of magnetic flux are 
arranged in a triangular Abrikosov lattice How- 
ever, thermal fluctuations destroy the flux lattice near 
the mean-field transition line and a flux liquid phase en- 
ters the phase diagram Q . As the temperature is reduced 
the vortex liquid undergoes a first order phase transition 
to what is commonly assumed to be the flux lattice state. 
This leads to a phase diagram as shown in Fig. |l|(a). The 
strong belief in first order melting of the Abrikosov flux 
latti ce r ests on the experimental evidence reviewed in 
Sec. 



[A 



Much of the analytical work on vortex lattice 
melting relies on the Lindemann criterion, which states 
that melting occurs if the mean fluctuation radius of a 
lattice point around its equilibrium position has reached 
a certain fraction (usually between 0.1 and 0.2) of the lat- 
tice constant Q. This criterion is not rigorous and does 
not provide a satisfying thermodynamic melting theory. 
The possibility of a first order phase transition due to 
decoupling of the different layers has also been investi- 
gated @. However, a decoupling transition is mostly 
expected to occur in addition to melting, and the lack 
of experimental evidence for two separate phase transi- 
tions has lead to a widespread belief that either there is 
no sharp decoupling transition or that it occurs simulta- 
neously with flux lattice melting. Our numerical results 
suggest a phase diagram which is fundamentally differ- 
ent from Fig. |l](a). It has a first order phase transition 



in excellent agreement with the first order transition line 
in YBa 2 Cu 3 7 (YBCO) in the B-T plane (see Sec. g). 
However, this transition is a decoupling transition and 
not associated with vortex lattice melting. There is only 
one phase in the phase diagram: the vortex liquid phase, 
and a vortex lattice exists only at zero temperature. 

Although there is striking experimental evidence for 
a first order phase transition in both YBCO and 
Bi2Sr2CaCu208 (BSCCO), there are certain features of 
the experimental data that are not explained by the stan- 
dard vortex lattice melting picture, most importantly the 
loss of first order behavior along the transition line at 
high (both for YBCO and BSCCO) and low (YBCO 
only) fields. Note that an end of the first order phase 
transition line at a critical end-point is not possible for a 
vortex lattice melting line, because the phase boundary 
separates phases of different symmetry. Our first order 
transition is not associated with any symmetry breaking. 
Thus the existence of a low field critical end-point should 
be expected and is directly observed in the simulation. 

In the framework of a vortex lattice melting picture 
the disappearance of the first order melting line can be 
explained by the presence of a tricritical point where the 
first order transition changes to a continuous one. Such 
behavior is commonly assumed to occur as an effect of 
sample disorder, which is to a certain degree present even 
in the best crystals. However, there is no wide consensus 
on the phase diagram in the presence of disorder. The 
three most important categories to distinguish are the 
disordered liquid, vortex glass and Bragg glass scenar- 
ios. For the first case there is no phase which is ther- 



1 



modynamically distinct from a vortex liquid and thus 
no thermodynamic phase transition. However, a fairly 
sharp crossover from fast to slow dynamics may occur 
within the vortex liquid state. The vortex glass scenario 
presented in detail in Ref. [|| relies on analogy to spin 
glass behavior. The vortex liquid is expected to freeze 
via a continuous transition to a vortex glass state char- 
acterized by short-range crystalline correlations but long- 
range phase correlations. Such a vortex glass phase would 
be truly superconducting with vanishing dc resistance. A 
popular recent theory predicts for weak disorder a first 
order transition to a Bragg glass state, which is charac- 
terized by slow, at most algebraic, decay of translational 
crystalline order |p|. 
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FIG. 1. Popular phase diagrams (a) in the clean case and 
(b) in the presence of disorder. Solid and dotted lines mark 
first order and continuous transitions respectively. 

With this inconclusive theoretical background, experi- 
mental and numerical evidence have had a major impact 
on the picture of the phase diagram of high tempera- 
ture superconductors (HTSC) including fluctuations and 
disorder. A popular phase diagram including disorder 
which accounts for many experimental features, notably 
the loss of first order behavior at high fields, is shown 
in Fig. 0(b). Our numerical results with disorder give a 
phase diagram with only one phase, a vortex liquid, just 
as in the clean limit. 



A. Experimental evidence 

This section attempts a review of experimental evi- 
dence on which both our and more conventional pic- 
tures of the phase diagram of layered superconductors 
are based. We discuss only evidence in YBCO, because 
this is the material to which our numerical model applies 
naturally. 



1. First order transition 

There is striking experimental evidence for a first order 
transition in YBCO. The earlier evidence for discontin- 
uous behavior suggesting a first order transition came 



from resistive measurements 0. A sharp drop in resis- 
tivity was found to occur at a temperature well below 
the H C 2 line. Later it was shown that these resistive 
drops coincide with a discontinuity in the magnetization, 
the first thermodynamic quantity found to be discontin- 
uous at the transition line The first measurements 
of the latent heat which unambiguously characterises a 
first order transition, were made by Schilling et al. |lC| ] 
in 1996. Since then a latent heat at the first order vor- 
tex transition has been observed in different crystals of 
YBa 2 Cu307_,5 with varying oxygen deficiencies 6 pi 12 
and for different orientations of the applied field . The 
B{T) dependence of the first order transition line obeys 
the standard continuum anisotropic scaling rules flifl un- 
der rotation of the applied field away from the c-axis [[l3) . 

The scaling behavior of the first order transition lines 
for samples with different oxygen deficiencies 8 and there- 
fore different mass anisotropies 7 is of some interest as 
it can be easily compared to predictions of different the- 
oretical models. In a range of samples the first order 
lines have been found to scale with 1/7 by Roulin et al. 
fL2fl , i.e. 7_B(T) collapses on one scaling curve. This is 
in disagreement with standard London-Lindemann type 
vortex lattice melting theory, which predicts the melting 
curve to scale inversly with the Ginsburg number Gi as 
1 /Gi cx 1 /"f 2 [|| . The 1 /7 scaling form is consistent with 
3D lowest Landau level (LLL) scaling and our numerical 
results. 



2. Loss of first order behavior 

The first order behavior at the vortex transition has 
been observed to vanish at an upper critical field B uc for 
samples which are not fully oxygenated Jl2[ | . B uc is found 
to increase with decreasing concentration of oxygen de- 
ficiencies S in YBa2Cu307_5. For 5 = a latent heat 
can be observed up to the highest experimentally inves- 
tigated fields of 16T The end of the first order line 
thus appears to be strongly correlated with the amount 
of point disorder in the form of oxygen vacancies present 
in the system. The upwards shift of B uc with increasing 
oxygen content fits in well with the theoretical phase di- 
agram in Fig. |l|(b), where it corresponds to an extension 
of the Bragg glass phase to higher fields with decreasing 
disorder. 

Another striking feature of the experimental first or- 
der transition line in YBCO is its termination at low 
magnetic fields, which has been consistently observed in 
all relevant calorimetric measurements |l^ [u|. The la- 
tent heat disappears for fields smaller than some lower 
critical field B\ c . The existence of a low field end-point 
is usually found "puzzling" Jl2| |. The variations of Bi c 
for specific heat measurements in different samples are 
large and qualitatively unexplained in the framework of 
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a vortex lattice melting picture. For B || c in near op- 
timally doped samples with a high level of oxygen defi- 
ciency S > 0.06, Bic is approximately 0.7T JlCj] . Mea- 
surements of Bic in different samples show that B[ c in- 
creases with oxygen content [ [L2| ], which suggests at first 
sight a correlation with twin density. The value of Bi c 
in fully oxygenated, twinned samples is of the order of 
several Tesla p| . The authors of Ref. |l2| discuss the 
origin of the end-point and the variation of Bi c in dif- 
ferent samples. The fact that detwinning does not no- 
ticeably change Bi c and the existence of Bi c in naturally 
untwinned samples, together with the reproducibility of 
Bic in different fully oxygenated samples, leads them to 
the conclusion that an intrinsic mechanism as a cause for 
the end point cannot be excluded. 

An important relation for our discussion of the value 
of Bic in different YBCO samples (see Sec. [II C ) is that 
an increase in oxygen content corresponds not only to an 
increase in natural twin density, but also to a systematic 
decrease in the anisotropy 7 in the samples used in Ref. 
^2| , ^5| . Our work suggests that this change in anisotropy 
rather than the presence of twins may cause the change in 
Bic. Our numerical work provides an explanation for the 
existence of Bi c as well as a qualitatively correct predic- 
tion of its rapid increase when 7 is decreased, such as can 
be achieved by increasing the oxygen content. Another 



noteworthy point which we shall discuss in Sec. [VB2 
is that the location of Bi c according to magnetization 
measurements is not always in agreement with the one 
measured in specific heat measurements. In a fully oxy- 



genated sample in Ref. 1 1 the latent heat vanishes at ca. 
6T while a magnetization discontinuity is still observed 
down to a field of 4T. 

Transport measurements reflect the loss of first order 
character of the transition for low as well as for high fields 
• The resistance only drops to zero, which would be 
the expected resistance for a weakly pinned lattice at the 
very low voltages used, for a limited range of magnetic 
fields. For high and low fields only a fractional drop is 
visible, which disappears completely somewhat below 2 
and above 7 Tesla. 

Below Bic and above B uc as well as in samples where 
no latent heat at all is observed, a "step" in the heat 
capacity C remains |l0|- |l^ , |l5|| . This behavior has been 
interpreted as evidence for a second order transition. The 
sharpness of this "step" , is not altogether convincing (see 
e.g. Fig. 3 in Ref. |b|]). However, the existence of a con- 
tinuous transition to a vortex glass state at high fields is 
expected for the theoretical phase diagram in Fig. |l|(b). 
According to the same phase diagram another line mark- 
ing a field driven phase transition line from Bragg glass 
to vortex glass is expected to emerge where first order 
melting turns continuous ||. The "fishtail" magnetiza- 
tion anomaly p^ , which correlates with the location of 
B uc [0 could be interpreted as evidence for such a tran- 
sition. A lack of sharpness of this feature makes it a 



candidate, however, for a crossover rather than a phase 
transition. There has also been evidence from resistive 
measurements for a field driven crossover line as an ex- 
tension of the first order transition in YBCO E§]. 



3. Evidence for a vortex lattice 

A vital ingredient of the vortex lattice melting scenario 
which this paper disputes is the existence of a vortex lat- 
tice. Experimentally a vortex lattice is indistinguishable 
from a liquid or glassy phase with short-range crystalline 
order on length scales large compared to the vortex sep- 
aration. Evidence for hexagonal coordination over large 
distances can be seen in YBCO for low fields in Bitter 
pattern decoration experiments fl9|| . At high fields this 
technique fails because the vortices are too close to be 
individually resolved. A powerful method to detect long- 
range vortex correlations is neutron scattering po| . The 
Bragg peaks observed in these experiments show that 
vortex positions are long-range correlated in all direc- 
tions. The correlation length along the field can be en- 
hanced by twin boundaries if the field is oriented along 
the c-axis. However, data from experiments with differ- 
ent orientation of the applied magnetic field show sim- 
ilar results, and thus indicate that the long correlation 
lengths along the field are independent of the presence of 
twin planes. The intrinsic crystalline in-plane correlation 
length is more difficult to deduce from neutron scatter- 
ing data, because twin boundaries and/or pinning to the 
underlying crystal determine preferred orientations and 
can thereby strongly enhance orientational order |2l| ] . 

Although neutron scattering experiments give evidence 
for long-range vortex correlations, some features of the 
data are unexpected in the framework of a vortex lattice 
melting picture. The observed diffraction patterns sug- 
gest the existence of a vortex lattice or a Bragg glass, 
which means that the melting transition is expected to 
be of first order. Such a first order melting transition 
should be visible as a discontinuous appearance of Bragg 
peak intensity, as the temperature is reduced. However, 
the peaks appear continuously, which leads the authors 
of Ref. [^0| to the conclusion that they must be dealing 
with second order vortex glass melting. 

An additional feature of the melting line defined by 
the onset of Bragg peaks (which is not mentioned in Ref. 
p(i|]) is that it lies in the B-T phase diagram distinctly 
below the line corresponding to thcrmodynamically mea- 
sured first order transition lines (l^Jf^l under the as- 
sumption of scaling with mass anisotropy like B cx I/7 
or B oc I/7 2 . This point will be investigated in more 
detail in Sec. IV A 2. For the alternative phase diagram 
presented in this paper, a continuous onset of Bragg peak 
intensity somewhat below the first order transition line 
is just the expected behavior. 
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B. Numerical simulations 



A large number of numerical Monte Carlo or Langevin 
dynamics simulations of different three dimensional mod- 
els have with few exceptions provided evidence of first 
order vortex lattice melting. The main disagreement be- 
tween different simulations of layered models concerns 
the question whether layer decoupling coincides with this 
melting transition. In this section we give an overview 
of different numerical results and point out what we con- 
sider to be the weaknesses of the respective models used. 

From the frustrated XY and Villain and lattice Lon- 
don models there is evidence for first order melting and 
distinct decoupling (22) as well as, at least in the thermo- 
dynamic limit, only one simultaneous first order melting 
and decoupling transition |2lJ. Very recent simulations 
of the uniformly frustrated 3D XY model show a vortex 
lattice melting transition as well as a second, possibly 
first order, phase transition within the liquid phase, at 
which the vortex line tension goes to zero (24). In these 
models vortices are confined to a lattice. In 3D, a lattice 
acts as a close grid of columnar pins with infinite pinning 
potential in the thermodynamic limit which may lead to 
spurious phase transitions. For the LLL, alteration of the 
phase diagram by the presence of a lattice pinning poten- 
tial which breaks translational and rotational invariance 
has been predicted from a theoretical analysis Pqj . 

There are also numerical models that avoid using a 
lattice. Numerical models relying on the 2D Bose gas 
analogy yield simultaneous melting and loss of phase co- 
herence along the c-axis |2(|. A different scenario has 
been seen in a simulation by Wilkin and Jensen [^7) , in 
which vortex pancakes in different layers are represented 
by particles with in-layer short-range repulsive and inter- 
layer attractive interactions. A first order transition as- 
sociated with decoupling of vortices and without melting 
character is observed. At a lower temperature a melting 
crossover without noticeable thermodynamic signature 
occurs. When sufficient point disorder is added, the de- 
coupling transition loses its first order character. While 
not being affected by pinning to a numerical lattice, the 
latter models may give unrealistic results because they 
allow variation of vortex position only, neglecting fluctu- 
ations of order parameter magnitude and in many cases 
having unrealistic short-range interactions psfl . Simula- 
tions using the Lawrence-Doniach (LD) model in the LLL 
limit, which allows for these fluctuations and which has 
long-range vortex interactions, show a single first order 
simultaneous melting and decoupling transition p9| , ^0[ . 
All of the simulations mentioned in this paragraph use 
periodic boundary conditions perpendicular to the field, 
which we believe can also lead to unphysical results (see 
Sec. g). 

In the following sections we introduce our numerical 
model (Sec. [Hj) and report results from our simulation. 



Comparisons to experimental data on YBCO are made 
in each section in the context of the relevant numerical 
results. Sec. Ill addresses the numerical phase diagram 



of layered systems in the clean limit, followed by an anal- 
ysis of order parameter correlations in space and time in 



Sec. IV. Numerical results in the presence of quenched 
random disorder are reported in Sec. [V} The paper closes 
with a discussion and a summary of our work. 



II. NUMERICAL MODEL 

Our simulation of a layered superconductor uses the 
Lawernce-Doniach (LD) model |n|, which consists of a 
stack of planes with Josephson coupling between neigh- 
boring layers. With the superconducting order parame- 
ter in the n th layer denoted as ip n , the Hamiltonian for 
the layered system in a magnetic field perpendicular to 
the layers is 



Hclean= ^ d J Sr (a\lp n \ 2 + ^p#n 



2m a b 



|(-iW-2eA)7/; n | 2 + J\ip n+1 -^„| ; 



where d is the layer periodicity, do is the layer thick- 
ness and B = V x A, which we shall take as constant 
and uniform. The same Hamiltonian can be read as the 
finite difference approximation to an anisotropic contin- 
uum model with tp(nd) = y d/doip n , (3 = fi^D xd/do and 
m c = fi 2 /2Jd 2 . In first approximation a(T) = a' (T — T c ) 
and (}2d{T) is constant; a' ,^d, J > 0. 

We simulate the LD model with N a b vortices per layer 
in N c layers. Along the c-axis we use periodic boundary 
conditions. In the afr-planes we chose a different, more 
unusual approximation to the thermodynamic limit of 
an infinite plane. The layers are taken to be of spherical 
geometry with a radial magnetic field. The reasons for 
our preference of this geometry to the more widely used 
geometry of a plane with periodic boundary conditions 
have been discussed in detail by Dodgson and Moore [B2| . 
The main advantage is that the spherical geometry guar- 
antees full rotational and translational symmetry, which 
periodic boundary conditions do not. One example where 
the spherical geometry captures the physics better than 
periodic boundary conditions - despite the topological 
defects imposed on the triangular lattice ground state on 



the sphere 



are particles interacting with the 1/r 



interaction. Here simulations on a sphere show already 
for moderate system sizes the genuine continuous tran- 
sition to the crystalline state ]33] |, while with periodic 
boundary conditions a spurious first order transition oc- 
curs even for very large system sizes. 

For each layer ip is expanded in eigenstates of the 
squared momentum operator (— iTiV — 2eA) 2 . We re- 
tain eigenstates belonging only to the lowest eigenvalue 
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2eBh, (the LLL approximation) which is a useful proce- 
dure over a large portion of the vortex liquid regime [j34| . 
Our numerical model is an extension to the model of a 
spherical thin film used in Ref . |32||35|,[36| . The magnetic 
potential is A = BR tan(#/2)</>, for which an orthonor- 
mal basis of the LLL eigenfunctions in each layer is given 
by 



An 



(AttR 2 ) 1 / 2 



l *sin m (0/2)cos JV ° b - m (0/2), (1) 



where A m = ((N ab + l)\/(N ab - to)!™!) 1 / 2 for < m < 
N a b. Note that we use units of length l m = (H/2eB) 1 / 2 , 
which fixes the sphere radius as R — (Nab/2) 1 / 2 . 

The order parameter in every layer is expanded in the 
above basis set as 



N 

■E' 

m=0 



(2) 



where Q = finksT/foodo) 1 / 4 - Orthonormality of the 
LLL eigenfunctions can now be used to express the 
Hamiltonian in terms of the LLL coefficients v n ^ m and 
only two parameters, U2T and to 



H 



clean 



N c / N ab ^ 2N ab 

X) ( "2T £ K, m | 2 + — Y, \Un,p\ 2 + 



n—1 \ m— 

N„. h 



J9=0 



\a2TV\ z2 \ V n+l,m ~ «n,m| 2 



(3) 



m=0 



with C/„ iP ({w„ ]m }) as defined in Sec. A 1, The 2D effective 
temperature and field parameter for each layer is olit = 
(doh/2e [32dB ksT) 1 ' 2 an , with an = ot(T) + eBh/m a b. 
The scaling parameter rj relates to the Josephson cou- 
pling constant J as rj = J/\au\- We can define an ef- 
fective mass to c via r\ — h 2 /2m c d 2 \aH |- Then ^/rj is the 
ratio of the 3D mean-field coherence length to the layer 
periodicity, ^/rj = £\\/d. Note that for a HTSC mate- 
rial the 2D parameters 02D and do are effective micro- 
scopic properties of the copper oxide layers and essen- 
tially unknown. However, they enter the simulation only 
via a.2T = (27rdo//32D ksTy^an where they can be re- 
placed by the layer periodicity d and the bulk (3 using 
the relation f3 = d/do x /?2D- 

The state in the LLL-LD model depends on two di- 
mensionless scaling parameters, «2T and 77. It is useful 
to have two scaling parameters which can be thought of 
as something physical, e.g. one characterising tempera- 
ture and the other coupling strength between layers. In 
this sense, «2T and ij are not very appropriate. Because 
c?o < d, the temperature parameter «2T goes to zero 
independently of B and T in the d — > limit of a contin- 
uous system. The coupling parameter 77 includes a factor 
|l/o!/f|, which means that it diverges at the mean-field 



H C 2. For these reasons we choose as effective tempera- 
ture and coupling strength two different parameters that 
depend on «2T and 77. 

For the temperature parameter in a layered model de- 
scribing a bulk sample, the 3D version of the LLL scaling 
variable ax stands out as an appropriate candidate. 
It is given by 



(XT - 



V2h 



3/2, 



2/3 



k B e 3 / 2 no 



1 

2 

K 7 



2/3 



dB c2 /dT\ Tc (T-T c )+B 



(BT) 2 / 3 



(4) 



which can be expressed in terms of «2T and r\ as \&t\ = 
?7(2a2T) 4 - At low temperatures |ar| 3 / 2 behaves as 1/T. 
A good measure of coupling strength is the product 
|a 2 T?7| which multiplies the coupling term in the Hamil- 
tonian. It is in SI units given by 



\Oi2Tm 



h 3 TT 



1/2 



8e 3 fc B AV n^d?/ 2 (BT) 1 / 2 ' 



(5) 



Other than the factor 1/ V BT, |a2T?7| contains only con- 
stants and therefore varies slowly for rather a wide range 
of an- It can thus be regarded as a material constant 
over considerable regions of the B-T plane. The limiting 
case |a2T??| — > describes the 2D system, while small 
|a2T J 7l means strongly layered characteristics. For con- 
stant ax, the limit |a2T 7 ?| ~~ * 00 is the continuum limit, 
in which all system properties depend on ar alone. Note 
that our model parameters depend on the bulk mate- 
rial parameters k, mass anisotropy 7, layer separation d, 
dB c2 /dT\ T = Tc , T c as well as B and T only. All of these 
are for HTSC more or less well known from experiment. 

Most of the experimental evidence to which we can 
compare our results is from samples which have at least 
weak disorder due to crystal defects or impurities. An 
idealized version of this kind of disorder is a random 
local variation of T c . To simulate the effects of such 
quenched random disorder, a random local potential O 
can be added to the quadratic energy term of every layer, 
to give a disorder contribution to the Hamiltonian of 



n dts =d J d 2 re(r)|?/;(r)| 2 , (6) 
where O(r) is real and Gaussian distributed with 



0(r) = 0, 



e(r)6(r') = AJ(r-r'). 



(7) 



Here S is the two dimensional Dirac delta function and 
A is the measure of the strength of the disorder. The 
random realizations of O are different in every layer, and 
Tidis is the sum of the single layer contributions. The dis- 
order contribution TLdis in terms of the LLL coefficients 
as used in our simulation is given in the Appendix, Sec. 
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A 2 . Unfortunately there is no experimental measure of 
A, and we can therefore not express the strength of dis- 
order in terms of measurable quantities. 

As in reference [[35), our simulation follows Langevin 
(model A) dynamics. We drive our system by the 
time dependent Ginsburg-Landau equation, discretized 
in time and expanded in the appropriate eigenfunctions: 

v n>m (t+At) - v ntm (t) = -Atr|^ + At^ m {t). (8) 

The complex random noise variables £„, m are drawn in- 
dependently from a Gaussian probability distribution, so 
that their magnitude has a variance a 2 / At = 1, where 
a 2 = 2 r fcflT, so that At is the only free parameter. We 
chose At = 0.15 (see Rcf. Q). 



III. NUMERICAL PHASE DIAGRAM 
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FIG. 2. Phase diagram (a) and magnetization discontinu- 
ity (b). Experimental data is taken from Ref. [^,[f0| (a) and 
Q (b). (For N c see section pIB| ) 



As an introduction to our numerical results we show in 
Fig. |^(a) our numerical phase diagram for the clean case 
in comparison to experimental results. The only phase 
present in this phase diagram is a vortex liquid. We see a 
first order transition line between two vortex liquid states 
with a critical end-point at low fields, which agrees well 
with the experimental YBCO "melting" line. The mag- 
netization jumps we observe are shown in Fig. |(b). They 
are in very good agreement with data for YBCO from 



Ref. gc}. The magnetization jumps in these experimen- 
tal measurements are very likely to be of thermodynamic 
origin, as they are found to be consistent (according to 
the Clausius Clapeyron relation) with the latent heat 
data measured in the same sample [[l0| . Figure ||(a) and 
(b) were obtained using standard YBCO values for the 
fitting parameters; viz for the Landau- Ginsburg param- 
eter k=60, the mass anisotropy 7=7.5, the slope of the 
mean-field transition line dB C 2/ '8T\t=t c = — 2T/K, the 
mean-field T C =92.5K and the layer separation d=11.4A. 



A. Continuum limit 



From LLL scaling we know that all thermody- 
namic properties depend on a 2 r alone in the 2D limit 
(|ck2T??| — * 0) and on olt alone in the continuum limit 
(|«2T??| — > oo). Figure ^ shows how the thermal average 
of a typical quantity of interest, here the Abrikosov ratio 
Pa = (|'0| 4 )/(IV'| 2 ) 2 i behaves in the intermediate regime 
of a positive finite coupling strength. 




|a 2T r||: 
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0C2T 

FIG 3. Abrikosov ratio (5a for a range of coupling 
strengths |Q2T??| plotted against the 3D and 2D effective tem- 
perature variables qt and ol^t- Solid lines are guide to the 
eye. 

In the high temperature regime 2D scaling applies, in 
the low temperature regime however 3D scaling becomes 
more appropriate. If we look at our model as a finite dif- 
ference approximation to the continuum case, this means 
that this approximation is for the same layer spacing 
d/(,\\ = 1/y^ better at low than at high 3D tempera- 
tures ax- This is a natural result as correlations along 
the field direction increase in the layered system as ax 
decreases. We are, however, not able to simulate sys- 
tem sizes that behave fully continuum-like at moderate 
temperatures. The numbers of layers used for the data 
in Fig. U are chosen to have the correlations along the c- 
axis distinctly smaller than the system size, which means 
for |a2T?7| = 12-8 even at the moderate temperature of 
a T = -6 an N c of 200. 
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B. First order transition 

Figure |](a) shows the phase diagram in terms of sim- 
ulation parameters, where data points mark the location 
of first order transition points as found upon heating 
and cooling the system. The logarithmic scale is cho- 
sen for even data distribution. As |ct2T?7| increases and 
the system approaches the continuum limit, the transi- 
tion line terminates at a critical end-point. Note that 
along the transition line ax is approximately constant, 
which means that the field and temperature dependence 
of the transition line behave as expected for a continuum 
model where ax is the only scaling parameter in the sys- 
tem. Another point of interest is the limit |a2T^| ~~ > 
of independent 2D layers. If we approximate the tran- 
sition value of ax as constant for different |a!2T?7|, this 
translates to a dependence of »2T on |a2T?7| along the 
transition line as ct2T oc — \ot2T'n\~ 1 ^ ■ As the coupling 
is reduced to zero, the 2D transition temperature di- 
verges, G-2T — * — oo. This would imply that there is no 
finite temperature phase transition in two dimensions, in 
agreement with 2D simulation results (3^j3(||35|] . Such 
behavior would be a natural for a transition of a pre- 
dominantly decoupling nature, which cannot occur in 
thin films. However, we find that direct extension of 
our numerical transition line to even lower couplings and 
thereby numerical observation of the 2D limiting behav- 
ior is impossible due to the finite size effects analyzed in 
Appendix O. 

Figures |J(b)-(d) show examples of the kind of measure- 
ment used to locate the first order transition. The system 
displays hysteresis upon heating and cooling. We mark 
the transition in the middle of the observed hysteresis 
loop. The hysteresis decreases with sweeping rate (typi- 
cally 10,000-20,000 time steps per data point), and for a 
few cases we have confirmed with equilibrium measure- 
ments that the equilibrium transition coincides roughly 
with the middle of the hysteresis loop. The coupling 
values for which we show hysteresis measurements corre- 
spond in decreasing order firstly to the critical end-point, 
secondly the nearest to the critical end-point where we 
have measured hysteresis, and thirdly an arbitrary, low 
coupling value. The discontinuity Ap in the order pa- 
rameter density p, given by p = (ax P/2nan) x (|V'| 2 )> 
is found to be more or less constant between |«2T?7| =0.1 
and 1.5. A decrease in Ap is observed below |at2T>7| =0.1. 
However, this decrease at low couplings is possibly due 
to finite size effects, which are in detail described in 
the Appendix ^|. The rapid decrease to zero between 
|ct2T??| =1.5 and |a2T??| =2.5 appears system size inde- 
pendent. 

The plots of the hysteresis in the degree of indepen- 
dence of neighboring layers, given by 

^ (|V(r)-^(r + dc)| 2 > 



reveal more about the nature of the transition and its dis- 
appearing. For low coupling, there is a large jump in T 
at the transition, which decreases throughout parameter 
space until it is very small just before the end-point. The 
decoupling character of the transition gradually decreases 
along the transition line until it disappears at the critical 
end-point. Further discussion predicting the existence 
and approximate location of the critical end-point from 
the decoupling character of the transition and more de- 
tailed numerical results concerning the critical end-point 
will be given in Sec. Ill C . The reader should note that 



in our simulation the onset of decoupling does not imply 
the onset of a superfluid density at the transition. 
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FIG. 4. Plot (a) shows the numerical phase diagram. First 
order transition points are plotted in the ar — \(X2TT]\ plane. 
Plots (b)-(d) show order parameter density p and the degree 
of layer independence T upon heating and cooling. Note the 
hysteresis in the system, the clear first order behavior for 
|a2T?7| = 2 (solid lines are linear fits of p above and be- 
low the transition) and the lack of first order behavior for 
|a2T??| = 2.5. N c varies between 8 and 80 for |a2T^| between 
0.02 and 2.5. For \a 2T r]\ = 2.5, p is offset by -0.05. 

The magnetization discontinuities in Fig. ||(b) are cal- 
culated from the discontinuities in p. The magnet iza- 
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tion in the LLL model is AttM — — (poeh/m a b)(\i^\ 2 ), 
which is in terms of our simulation parameters AirM = 
tt(B — B c2 (T))p/ctTK 2 , where B is the applied magnetic 
field. Thus we can work out the magnetization disconti- 
nuity from the discontinuity in p taken from the two lin- 
ear extrapolations at the transition. The data points in 
Fig. |(b) represent \a 2T v\ =1, 1-5, 2 and 2.5 for N c =50, 
60, 60/80, 80. For these four transition points we have an 
average value of cut — —7.72, which yields the transition 
line in Fig. p|. 

The magnetization discontinuity is thermodynamically 
linked to the entropy change at the transition by the 
Clausius-Clapeyron relation 



AS = -AM 



dHpoT 
dT 



(9) 



where dHpor / dT is the slope of the first order transi- 
tion line. We cannot measure AS directly. However, 
the good agreement in AM and the first order transi- 
tion line between simulation and experiment on the one 
hand and the consistency of the experimental AS as cal- 
culated from magnetization and latent measurements for 
the samples in Ref. on the other hand imply agree- 
ment for AS between simulation and experiment. Note 
also the clear change in slope of the linear fits to p at the 
transition. Locally cxt is linear in T, and via Maxwell's 
relation 



dS 
dH 



dM 
~dT 



H 



the sudden change in slope of p at the transition implies 
a change in slope of the entropy and a step-like feature in 
the specific heat C = T(dS/dT) B , with a lower value on 
the low temperature side of the transition. This is con- 
sistent with experiment. The relative change in slope we 
find from the fits in Fig. [|is 8%, while the equivalent ex- 
perimental change in the heat capacity as taken from Fig. 
3 of reference [|l3| is of the order of 5% (for a derivation 
of this value see Sec pIFj ). The fact that our simulation 
gives evidence for a step in the specific heat is consistent 
with results from a theoretical analysis showing that the 
step seen in experiment can be accounted for by thermal 
fluctuations within the LLL approximation feTf. 



C. The critical end-point 

Until now we have as evidence for a critical end-point 
only the fact that the hysteresis along the transition line 
eventually becomes unobservable. To more firmly estab- 
lish its existence, we shall consider how the nature of the 
transition may lead to a critical end-point. 



1. Why does the first order transition disappear? 

A well known first order phase transition which ends 
at a critical end-point is that of the ordinary liquid-gas 
transition. Here the phase transition separates a liquid 
state with small interparticle separation di which takes 
advantage of the attractive interparticle energy which ex- 
ist at distances d m i n , so di » d m i n , from a gaseous state 
with large interparticle separations d g which is favored 
by a high entropy. If the density is increased to the point 
where it reduces d g to be of order d m im the transition 
line ends. We believe that in our case the entropy advan- 
tage of the high temperature phase arises when the order 
parameter values in adjacent layers are uncorrclatcd i.e 
when the layers are decoupled. The ratio of the mean- 
field coherence length perpendicular to the layers to the 
layer distance, £,\\/d = y/rj, increases along the transition 
line with increasing coupling parameter |«2T??|- Because 
£|| defines the minimal extent of order parameter cor- 
relations, a high-entropy state with decoupled layers is 
not possible if £n > d. And indeed, we will in the next 
section locate the critical end-point where the transition 
disappears at |a2T??| = 2.55 and ay = - 7.75, which corre- 
sponds to y/rj — 1.06. Very near the zero-field transition 
temperature T c , where £| S> d, the system can be ex- 
pected to behave like a continuum and thus a decoupling 
line cannot be expected to reach T c . 



2. Divergence of length scales 

Near a critical end-point we do not only expect all dis- 
continuities to disappear, but we also expect there to be a 
divergence of the length scale of fluctuations in the order 
parameter density of the system. We therefore looked at 
the density-density correlations of the order parameter, 



C d (r',t') = 

Mr, t)| 2 |V(r + r',t + t')\ 2 ) - (|<Mr)| 2 ) (|^(r + r')| 



(10) 

2 > 



M 2 ) 2 

in the case where r' is a vector parallel to the c-axis and 
t' = 0. This correlator is expressed in terms of thermal 



averages of the LLL coefficients in the appendix, Sec. B 2 . 
Plots of these correlations near the critical end-point can 
be seen in Fig. ||. 

There is evidence of two length scales in the vicinity 
of the end-point. The short distance decay of the cor- 
relation function is dominated by the positional correla- 
tions of the vortices in the different layers. This length 
scale is mostly determined by ax and changes slowly 
in the vicinity of the critical end-point. The diverging 
length scale is thus not that of positional correlations of 
the vortices, but instead is associated with local density 
fluctuations. This is not surprising if the analogy of an 
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ordinary liquid-gas transition is considered, where crys- 
talline correlations on a microscopic scale correspond to 
positional vortex correlations. A noteworthy feature of 
the short distance decay of vortex correlations is that the 
difference between curves at olt = —7.6 and ar = —7.8 
decreases as the coupling parameter |a2T»?| is raised past 
the end-point, as one would expect for a disappearing 
discontinuity. 
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FIG. 5. Density-density correlations along the c-axis near 
the end the first order transition line and the critical end-point 
(see inset). Note the decreasing difference in correlations be- 
tween ot = —7.6 and «t = —7.8 as the transition disappears. 
For |q2t??| = 2.5 and ar = —7.8 we see evidence for a long 
length scale associated with fluctuations in the average order 
parameter density. System sizes are N a b=72, 220 < N c < 260 
for a T = -7.6 and N c = 270 for a T = -7.8. 



Density correlations diverge at the liquid-gas critical 
end-point on a mesoscopic scale. Analogous correlations 
on a larger scale appear in the numerical data between 
ax = —7.6 and cxt = —7.8 as laurel i s increased to its 
critical end-point value. Figure [| shows for |«2T'7| =2.5 
and o>t = —7.8 evidence of a second, much longer length 
scale governing the decay of the correlation function at 
large distances. This length scale is associated with the 
density fluctuations at the critical end-point and only 
becomes visible once it is larger than that of the vortex 
correlations. Due to the small amplitude of these density 
fluctuations very long simulation times are needed to see 
the correlations within the statistical noise. 



D. Anisotropic scaling and the value of Bi c 

The location of the numerical first order transition is 
to first approximation a line of constant ax and thus in 
agreement with 3D LLL scaling. This is, although sur- 
prising for very low couplings, not unnatural in YBCO, 
where even for the highest fields and lowest couplings, e.g 
for a transition at B=20T, T=65K, \cxtt]\ is of the order 
1/2 and the correlation length along the c-axis above the 



transition (see Sec. IV A I) of the order of 10 layers. The 
continuum approach can thus be expected to work fairly 
well for YBCO. 



1. Scaling with anisotropy 7 

According to Eq. || the BIT) dependence of a line of 
constant ay ban be approximated as 



B 



K 2 7 



Thus we expect our transition line in samples of differ- 
ent anisotropies to scale as B oc 1/7, which is the ex- 
perimentally observed scaling p2[ . This form of scaling 
disagrees with a London-Lindemann type melting theory. 
The variation of the location of the low field critical end- 
point with anisotropy 7 is discussed at the end of this 
section in the context of a general analysis of variations 
in B ic . 



2. Variation with the angle of the applied magnetic field 

We cannot in our simulation change the orientation of 
the applied field. We can, however, using the known 
properties of the numerical transition discuss the ex- 
pected behavior upon such a change in field orientation. 
We believe that the first order transition is predomi- 
nantly of a decoupling nature. However, it is important 
to keep in mind that the coupling between layers in our 
simulation is not magnetic coupling of vortex segments, 
but Josephson coupling of the order parameter. This 
type of coupling is independent of the orientation of the 
flux lines with respect to the c-axis. Under angle rotation 
between = (B\\c) and 9 = ir/2 (Sic), the transition 
line should scale according to anisotropic continuum scal- 
ing Q with 6 as 

B(8,T)= 7e B(0,T), 

where 70 is given by jg = (cos 2 6* + sin 2 6> / 7 2 ) -1 ^ 2 . This 
form of scaling has been observed for the first order tran- 
sition line in YBCO @. 

Our argument linking the nature of the transition to 
the location of the end-point applies for any orientation 
of the magnetic field. The location of the end-point is 
for all orientations given by £ c ~ d, where £ c is the co- 
herence length along the c-axis. For B\\c this condition 
is equivalent to £|| w d. Using £|| = h/(2m c \aH\) 1 ^ 2 , 
this condition can be transformed to a simple B(T) de- 
pendence which scales like the transition line itself. The 
end-point where both lines cross therefore equally just 
shifts to a higher field as Bi c (0) = jeBi c (G). This form 
of scaling has been experimentally observed for the end- 
point in YBCO Q. 
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3. Variations in Bi c 



The location of the numerical critical end-point agrees 
well with the experimental data of Fig. However, as 



already mentioned in Sec. I A, the experimental value 
of Bi c varies widely between different samples. Approx- 
imating k and T as constant at the critical end-point, 
we obtain from Eq. [| the scaling relation B\ c oc I/7 4 . 
The K dependence is Bi c oc 1/k 2 . These two scaling rela- 
tions show that Bic depends very sensitively on material 
parameters. 

We can use the fit in Fig. |with 7 = 7.5 and B ic = 0.7T 
as reference point to compare the location of the critical 
end-point in our simulation to experimental values of Bi c . 
For samples with measured anisotropics of 7=7.8 (Ref. 
fL3f), 7=7.0, 5.9, 5.3 (Ref. jl5)) specific heat peaks have 
been measured down to B Zc =0.7T (Ref. (l3)), B ic =3T, 
4.5T, 6T (Ref. @). The rapid increase of B lc with de- 
creasing anisotropy is in qualitative agreement with a 
scaling law Bi c oc I/7 4 . 

Quantitative analysis however yields only poor agree- 
ment. The predictions using the above reference point 
and scaling law deviate from the experimental value by 
-15% (Ref. ||) and -70%, -60%, -50% (Ref. 
There are many possible reasons other than variations in 
k for this quantitative disagreement. Firstly, the finite 
width of the transition due to sample inhomogeneities 
may lead to a spreading out of the specific heat peak, 
which can make it undetectable for the lowest fields above 
Bic and lead to overestimates of Bi c . Secondly any as- 
pects of the physical coupling which are not represented 
in our model could lead to corrections in the effective 
|c*2T?7| and should thus be included for an accurate de- 
scription. 

A third, very important point is that near T c critical 
fluctuation effects arising from the zero field transition 
are not negligible and especially affect the divergence of 
Such effects extend to fields of the order of Gi x H C 2 (0) 
|§, in YBCO - IT. Up to these fields, the LLL ap- 
proximation is invalid because higher Landau levels are 
needed to allow for critical fluctuations. Thus the end- 
point lies in a region where the LLL approximation is 
inadequate, and we cannot expect our estimates of the 
position of the critical end-point to be quantitatively ac- 
curate. However, we can expect that such fluctuations 
well make Bi c smaller than estimated from the LLL ap- 
proximation. Critical fluctuations cause £|| to diverge 
faster than in the LLL approximation, so it reaches a 
value O(d) at temperatures closer to T c than in the LLL 
approximation. This in turn will reduce the value of Bi c . 

In addition, despite extensive finite size effect analysis 
(see appendix ^) we can never fully exclude the possi- 
bility that the location of the end-point would shift to 
lower fields if we used much larger systems. The exis- 
tence of an end-point is fairly securely established by our 



physical argument. We have also made sure that the the 
ratio of correlation lengths to the linear system dimen- 
sions, l a b I Lab and l c /L c , decrease rather than increase 
as we pass the critical end-point with increasing |a2T??|- 
However, the general tendency of small system sizes to 
decrease discontinuities at the transition could have lead 
to an overall underestimate of the jumps due to finite size 
effects and thus to an overestimate of Bi c . 

The presence of sample disorder increases the value of 
Bic, an effect which we report for random disorder in 
our simulation in Sec. Real YBCO samples exhibit 
in general not random disorder but clusters of oxygen 
vacancies and large scale sample inhomogeneities. Oxy- 
gen clusters pin at low fields a considerable fraction of 
the vortex matter, so that the field can be divided as 
B = B p i nne d + B f r ee- The fraction of pinned vortices will 
be increased by the presence of sample inhomogeneities 
which provide regions at an effectively lower olt that ex- 
hibit stronger pinning. Such a decrease in the effective 
field of the free vortices (Bf ree < B) could considerably 
increase the observed Bi c . 

Considering that the magnetic field at the end-point 
depends so sensitively on the model parameters, critical 
fluctuations and disorder, the good quantitative agree- 
ment with the data of Ref. Jn],|l3| seems more than can 
be reasonably expected from our model, and maybe is 
indeed a product of chance in which Bi c is in our simula- 
tion increased by inaccuracies of our model and/or finite 
size effects by the same amount that it is increased by 
disorder effects in the samples in Ref. ]To|,^3[. 



E. Beyond the critical end-point 

It is often supposed that the first order transition in 
YBCO changes to second order below the end-point, 
where no latent heat is visible but a "step" in the heat 
capacity C remains [[l3],[l2| . We believe however that this 
"step" can be identified with the onset of a small rounded 
peak in the superconducting specific heat C s = C — C„, 
(n for normal state), which is known to arise from ther- 
mal fluctuations |39j| . In this section we shall examine 
experimental specific heat data and numerical data in 
the light of this possibility. 

The specific heat peak due to thermal fluctuations has 
been observed for example in niobium by Farrant and 
Gough [^(J, where observations are in good agreement 
with theoretical predictions |$9|. We find that the loca- 
tion and the height of the peak as well as the length 
of the rise (or width of the "step") in C from the 
low temperature value C Stm f (mf for mean-field) to its 
maximum agree well for the niobium and YBCO mea- 
surements taken from Ref. ^] and p3| . We shall now 
explain in detail how the data from both references can 
be compared and give respectively values for the loca- 
tion, width and height of the specific heat peak or "step" . 
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For a clearer picture of how the experimentally measured 
specific heat splits into normal and superconducting con- 
tribution as well as how it compares to the mean-field 
contribution, a schematic plot is given in Fig. 0. 



C=C s +C n 






""\"""" 

C n 



FIG. 6. Schematic view of a typical experimental specific 
heat curve without first order peak (thick line) with normal 
state contribution (dashed line) and mean-field value (thin 
line) . 

Farrant and Gough give in Ref. the superconduct- 
ing specific heat data already in terms of LLL scaling pa- 
rameters. The plotted quantity is C s /C Sjm f . The data 
in Fig. 3 from Ref. |l3) which shows the specific heat in 
YBCO is given as C minus C(B — 0). The latter near 
the "step" is approximately equal to the low temperature 
value C s>m f + C n , because for B — ► 0, qt-* — oo. This 
approximate equality is also visible in Fig. 1 of the same 
Ref. jD§. The plotted quantity C-C(B = 0), is therefore 
approximately equal to C s — C s>m f. 

The peak in C s in niobium obeys LLL scaling and 
is found at cut ~ —7. The specific heat maximum in 
YBCO occurs at temperatures just above the extrapo- 
lation of the first order transition line, which is located 
at a>T ~ —7.8. For the example curves for B — 0.25T 
and B = 0.5 in Fig. 3 from Ref. the center of the 
broad specific heat peak is at T w 91. 4K and T sa 90. 7K 
respectively, which both translate with the same previ- 
ously used YBCO material constants to ar — —7.2. This 
is in very good agreement with niobium. 

The width of the specific heat rise in niobium A«t ~ 2. 
For B = 0.25T, no sharp step feature is visible in the 
YBCO data. The specific heat rise from the low temper- 
ature value to the maximum takes place in the tempera- 
ture region 91 — 91. 4K, which corresponds to Aar = 2.8, 
broader than in niobium. For B — 0.5T one might sus- 
pect a step- like feature located 90 — 90. 5K. This width 
corresponds to Acvt = 2.3, a value in agreement with the 
niobium data. 

For niobium C s is at its maximum 5% larger than 
C Sy mf. In YBCO we have to divide the plotted data 
by C Sim f to compare with this value. C s>m f is roughly 
given by the step in C at the zero field transition, 
which we take from Figure 1 of Ref. |b|] as C Sim f w 
60 mJ/mole K 2 . For all fields the specific heat "step" 
in YBCO is of the order of 1.5 mJ/mole K 2 , which gives 
(C s — C s ^ m f)/C Stm f Ri 2.5%, which we consider as reason- 
able agreement with in niobium. Exact agreement cannot 



be expected for the following reasons. The value of C s . m / 
we used for YBCO is a rather crude approximation. For 
niobium C S)m f is also somewhat uncertain, because it 
depends on the choice of C n (see Fig. which is not 
directly measurable but extrapolated from fits to higher 
temperature data. Also the data in YBCO at fields be- 
low Bic are in a region where the LLL approximation is 
inadequate. 




FIG. 7. Plots of the order parameter density p upon heat- 
ing and cooling for «2T??| — 3. Solid lines are linear fits for 
the regions above and below the extrapolation of the first or- 
der transition line, extrapolated as dotted lines. N a t = 72, 
N c = 80. 

If the step feature is mainly due to the increase of LLL 
fluctuations, it should be observable in our simulation. A 
"step" in the superconducting specific heat corresponds 
to a change in the slope of the magnetization, in simula- 
tion terms a change of dp/dotT- Figure shows p upon 
cooling and heating for |ct!2T?7| — 3, well beyond the crit- 
ical end-point. The data for or < —7.9, i.e. below the 
region where the first order transition takes place in more 
layered samples, is very strongly affected by hysteresis. 
The sweep rate being of the same order as usual, this is a 
sign of very long fluctuation time scales. We shall come 
back to this point in Sec. 



[V. Linear fits to the data for 



both cooling and heating below and above the extrapo- 
lated transition line, ax < —7.9 and «t > —7.8, give a 
change in slope of 8%, larger than in the experimental 
data from both niobium and YBCO. 

In summary we have seen that the specific heat "step" 
in YBCO at different fields has approximately the same 
amplitude as well as width and position when expressed 
in terms of ax, i-e. the "step" feature obeys LLL scal- 
ing (LLL scaling has for the steps associated with the 
first order transition previously been established in Ref. 
p7{). The corresponding data from our numerical simu- 
lation are not equilibrated and therefore not very accu- 
rate, but consistent with a similar "step" feature. The 
semi-quantitative agreement between YBCO and nio- 
bium strongly suggests that we are dealing with the same 
phenomenon and therefore that there really is no sharp 
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specific heat step in YBCO beyond the critical end-point. 



IV. ORDER PARAMETER CORRELATIONS 

The existence of the critical end-point implies that no 
symmetries are broken at the transition, which means 
it cannot be a liquid-crystal transition. Investigation of 
the nature of order parameter correlations described in 
this section confirms that the state below the transition 
is a vortex liquid. However, neutron diffraction patterns 
corresponding to a triangular lattice and an electrical re- 
sistance close to zero below the first order transition line 
suggest very long-range vortex correlations. We shall ar- 
gue that the extremely fast growth of correlation length 
scales and relaxation time scales below the first order 
transition, which have been theoretically predicted in the 
low temperature regime f4lf| as well as observed in our 
simulation, can account for these effects. 



A. Static order parameter correlations 

We have measured various different kinds of equilib- 
rium order parameter correlation functions, including the 



structure factor as previously defined in Ref. |32| and 
given by the normalized density-density correlator in Eq. 
|ll] (below) for Ar = along the c-axis. The peaks in 
the structure factor, which occur at the reciprocal lattice 
vectors of the triangular lattice, reflect the crystalline 
correlations within each layer. To examine order parame- 
ter correlations along the c-axis, we have measured three 
different correlators. Firstly, we measured the density- 
density correlations along a line parallel to the magnetic 
field, Cd(Ar || c), the same correlator that gave in Sec. 



Ill C evidence of a diverging length scale at the critical 
end-point. Secondly, we measured the decay of the ab- 
Fourier transformed density-density correlations, which 
are for Ar || c = the structure factor of the system, 
depending on Ar || c. Just as for the structure factor J3^] 
we normalize by the high temperature limit, so that the 
relevant quantity is 



C d (k\\ab,Ar\\c)/ lim C d (k\\ab,0) , 

OLlT — 



(11) 



with 



C d (k||a6,Ar||c 
1 



, |(i|22 «kf(k, r)H 2 (-k,r + Ar))- 

(H 2 (k, r))(|Vf(-k,r + Ar». 

The definition of this correlator in its generalized, time 
dependent form in terms of the LLL coefficients and the 
high temperature limit are given in Sec. 



B2 



Eq. B2 



and 

B3|. We find that this correlator decays most slowly with 



Ar || c for k a b ~ G, the first reciprocal lattice vector. In 
addition to Cd we also measured the phase correlator C p 
defined as 



C p (Ar||c) = a T 



/3(V>*(r)y>(r+Ar)) 

27TCt£f 



(12) 



which is expressed in terms of the LLL coefficients in Sec. 
B2. The prefactor is chosen such that that C p (0) is the 
order parameter density p. In Fig. || we show examples 
of static density and phase correlations above and be- 
low the transition for |ce2T^| = \ which corresponds to a 
transition temperature of 83 K in YBCO. 
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FIG. 8. Static order parameter correlations above and be- 
low the first order phase transition at |a2T??|=l. (a) Structure 
factor (with Lorentzian fits), (b) phase correlations along the 
c-axis. Plots (c) and (d) show density-density correlations 
along the c-axis, where the correlator in (c) is in real space 
for points with the same coordinates in the a&-plane, and the 
correlator in (d) is in Fourier space near the first reciprocal 
lattice vector of the triangular lattice, G. Note that in plots 
(b),(c),(d) the y-axes are such that parallel linear fits corre- 
spond to the same length scale. The growth of length scales 
is slowest for the phase correlations in (b) and fastest in (d). 
The system size is N ab =72, iV c =270. 

Figure||(a) shows examples of the structure factor with 
Lorentzian fits. The inverse width at half maximum, 
(5" 1 , of the peak near the first reciprocal lattice vec- 
tor is proportional to the crystalline correlation length 
l a b within one layer [ p2[ , S^ 1 = l ao /2l m . No qualita- 
tive change in behavior is visible across the transition, 
just sharpening of the peak, which corresponds to an in- 
crease of l ao . For the above fits l a b/lm =2.66, 3.14, 3.88, 
4.94, 5.46. This is shorter than the radius of the sphere 
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R/lm — y/N/2 — 6 and shorter than the average dis- 
tance between the 12 disclinations imposed by the spher- 
ical geometry, y/ 4irR 2 /l2 « R. This means that finite 
size effects on this data are not expected to be too large. 

Figures ||(b)-(d) show examples of measurements of 
the three different correlators along the c-axis. Figure 
|8|(b) shows the phase correlator, Fig. |§|(c) the real space 
density-density correlator and Fig. |£](d) the a6-Fourier 
transformed density-density correlator near the first re- 
ciprocal lattice vector G. In all cases we see an exponen- 
tial decay of the correlation function with a finite length 
scale l c for density correlations below as well as above 
the transition. Only for a liquid phase would these cor- 
relation functions all have an exponential decay. 

All three correlators decay over half the system size 
ri < N c /2, and then rise to C{N c d) = C(0) due to 
the periodic boundary conditions. To extract decay 
length scales we fit to the first portion, approximately 
Ar < (N c d/6), of the decay, as indicated by linear fits. 
This region of decay is least affected by finite size effects, 
statistical noise and errors due to incomplete equilibra- 
tion (see Sec. ^J). The length scales extracted from the 
linear fits are of the same order for all three correlators. 
Just above the transition the different l c are almost equal. 
Below the transition the Fourier transformed density- 
density correlations clearly have the largest length scale 
with l c ,d/(,\\ = 58 (d for density) at ax — —8.1 and the 
phase correlator the smallest with l c ,p/£\\ — 42 (p for 
phase) at the same ax- 

At temperatures well above the transition we find 
the opposite behavior. The density-density correlation 
length is at temperatures ax ~ —4 little more than 
half of the phase correlation length, which agrees with 
the simple high temperature expectation (ip(r)*ip(r')) 2 ~ 
(| , 0(r)| 2 |'!/'(r')| 2 ) and therefore for exponential decay 
l cp rj 2l c _ d- However, because we find that below 
the transition l c ^d is the longest and therefore dominant 
length scale, we shall in the following refer to the decay 
length scale of the correlator Cd{k a b = G) as l c . 



1. Temperature dependence of correlation lengths 

In Fig. U we plot on the left the ax dependence of 
the inverse width at half maximum of the first peak in 
the structure factor, (5" 1 , which is proportional to the 
length scale of in-plane crystalline order, l a b/lm = 2<5 _1 , 
and on the right the length scale l c as obtained from lin- 
ear fits to the decay of density-density correlations. The 
length scales have a discontinuity at the transition, which 
is found to grow with distance from the end-point, as 
one would naturally expect. This discontinuity is clearly 
visible for |ev2T?7l=0.05. At large couplings, finite size ef- 
fects spread out the discontinuity. The data for |a2T^|=l 
shows a rapid growth of the length scales at and below 
the transition. 



If the phase coherence of the Abrikosov state is exam- 
ined in the presence of thermal fluctuations, one finds 
that in and below three dimensions thermal fluctuations 
destroy phase coherence at any finite temperature [ fL2| . 
Under the assumption that there is only one relevant 
length scale describing phase order and crystalline order 
in the system, perturbative studies for the continuum 
low temperature regime predict an exponential growth 
of length scales with |cvt| 3 / 2 as l c oc exp(^4|o!T| 3 ^ 2 ) and 
l ab cx exp(0.5A|a T | 3 / 2 ) glj. The authors of Ref. @ 
estimate that A may be given by its upper limit value 
A = 0.53. The slope of such growth behavior with 
A = 0.53 is given by the solid lines in Fig. [| The growth 
rate in our simulation data is at the lowest temperatures 
of the same order as the theoretically predicted upper 
limit. However, the analytical result is from an expan- 
sion around zero temperature in the continuum limit. 
This regime is for numerical reasons described in Sec. ^ 
not accessible to our simulation. Therefore perfect agree- 
ment of our simulation results with this form cannot be 
expected. 
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FIG. 9. Logarithm of correlation lengths plotted against 
|o!t| 3 ^ 2 . Solid lines represent the growth rate for the low tem- 
perature regime as predicted from perturbative expansions 
around zero temperature Dotted and dot-dashed lines 

mark the width of the first order phase transition from the 
magnetization discontinuity for |«2T?7|=0.05 and ja2T?/| = l re- 
spectively. Note the fast growth of length scales for |a2TT?|=l 
below the phase transition. Note also that the data below the 
transition for |a2T^j=0.05 is strongly affected by finite size 
effects (l a b ~ 10Z m > R = 6/ m ). 

The data in Fig. | for |a2T7?|=3.2, which is beyond 
the critical end-point and has got no phase transition, 
suggests that a faster growth of length scales sets in ap- 
proximately where the phase transition is located in more 
layered samples. The reason why we have obtained no 
more data at lower temperatures to confirm this tendency 
is just because of the fact that length and time scales 



(see Sec. EVB) grow so fast that we found it impossible 
to equilibrate a sufficiently large system at even lower 
temperatures. (Remember that the same system size in 
units of £|| corresponds to more layers for larger |cv2T?7|-) 
However, the onset of fast growth of length scales can 
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be expected from the fact that in the low temperature 
regime, where length scales are so long that the contin- 
uum approximation works well even for strongly layered 
materials, the data for different |a2T?7| must collapse onto 
a single curve and depend only on olt- At the transition, 
length scales in the more continuous samples are in the 
appropriate scaling units distinctly shorter than for low 
|a 2 T??|- For the length scales in these systems to approach 
the lengths scales of the more layered systems, rather fast 
growth is necessary not far below the transition. 



2. The appearance of Bragg peaks 

Just below the phase transition the correlation length 
scales in the simulation are for the range of coupling pa- 
rameters that correspond to YBCO not comparable to 
the much larger length scales needed to give a signal in 
neutron diffraction experiments. Although the structure 
factor in a liquid is rotationally symmetric, coupling with 
the underlying lattice or preferred orientations given by 
twin planes may for long length scales lead to the ap- 
pearance of Bragg- like peaks Q| . 

Our simulation suggests that the vortex liquid is not 
far below the first order transition correlated and effec- 
tively crystalline over length scales comparable to the 
system size or a "Larkin"-like length scale (dependent 
on the amount of disorder present). We can extrapolate 
the growth in length scales in Fig. |^ below the transition 
assuming the exponential growth rate estimated in Ref. 
J4l"| and indicated by the solid lines. For a decrease of 
AaT « 1.2 for example, which corresponds in YBCO to 
cooling by only IK or (1/4)K below the transition at 5T 
or 0.7T respectively, we obtain an increase in l a b by a 
factor of 4 and an increase in l c by a factor of 16. At 
ax ~ —10.5 the crystalline correlation length l a b has ac- 
cording to the same estimate reached 30 lattice spacings 
and l c the order of 10,000 £||, which is for magnetic fields 
of 5T or 0.7T of the order of 5,000 and 10,000 lay- 
ers respectively. Correlation lengths of this order can be 
expected to lead in real samples to observable neutron 
diffraction peaks. 

Figure ^ shows scaling plots of different experimen- 
tal phase diagrams in different YBCO samples. We ex- 
pect according to LLL scaling that the plots of jB(T) 
from samples with a different mass anisotropy but oth- 
erwise identical material parameters collapse (see Sec. 
HIBl). Roulin et 



Jl2| report that the first order tran- 
sition lines in different clean samples, represented by a 
solid line in Fig. fo, as well as the the line of steepest 
slope at the specific heat step in disordered samples [[l5| , 
represented by a dashed line in Fig. |l0L collapse in this 
way. The solid and the dashed line correspond to the 
specific heat peaks in three different samples with differ- 
ent anisotropies each. The Schilling data (scaled using 



7 = 7.8 |r|) also agrees reasonably well with this scaling 
form. 

According to the vortex lattice melting picture, the 
points which mark the onset of neutron diffraction peaks 
in Fig. [n] (scaled using 7 = 4.3 p(}]), should coincide 
with the first order transition line. Comparison however 
shows that is roughly a factor 2 below the first order 
transition lines. (In a scaling plot of 7 2 -B, which should 
yield a collapse assuming London-Lindcmann-type melt- 
ing, the discrepancy between the Aegerter and Schilling 
data would be even larger, roughly a factor 4.) Very sub- 
stantial differences in other material parameters would 
be necessary to make up for such a large deviation. In 
our picture the relative position of the two features is 
natural, because the onset of neutron diffraction peaks is 
attributed to a crossover below the first order transition, 
when the neutron diffraction experiment becomes sensi- 
tive to the exponentially fast diverging length scales. 
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FIG. 10. Comparison of phase diagrams obtained on dif- 
ferent samples from thermodynamic and neutron diffraction 
measurements. All data is for B || c. Empty circles and the 
solid line (extrapolated to T c beyond Bi c ) correspond to spe- 
cific heat peaks. The dashed line represents the point of steep- 
est slope of the specific heat "step" attributed to second order 
melting in samples which do not show a specific heat peak. 
The data for the solid line, empty circles, the dashed line and 
filled circles is taken from Ref. |12| , |l3|jl5[|20| respectively. 

However, one might say that comparing first order 
"melting" and onset of neutron diffraction is not a com- 
parison of like with like, because the sample in the neu- 
tron diffraction experiment is far too dirty to exhibit a 
first order phase transition. The authors of Ref. |2(| es- 
timate that their sample is comparable to samples from 
Ref. [jl5| and attribute the onset of neutron diffraction 
peaks to second order freezing to a vortex glass. In this 
case they should coincide with the line marking position 
of the steepest slope of the specific heat steps in such 
samples, which has been interpreted as a second order 
vortex glass transition However, the plots in Fig. [Toj 
show that this is not the case. 
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Our scaling approach is based on the approximation 
c*t ~ l/(7-S) 2 ^ 3 - Assuming that the first order transi- 
tion line found by Roulin et al. [ fl2[ occurs as in our simu- 
lation at ax ~ —7.5, this means that the dashed line and 



the onset of neutron diffraction in Fig. 10, for which 7B 
is reduced by a factor of about 3/4 and 1/2 respectively, 
correspond to out ~ —9, and o>t ~ —10.5 respectively. 
The olt value for the onset of neutron diffraction peaks 
is in excellent agreement with our estimates as to where 
l c « 10,000£|| and l a b ~ 30 lattice spacings in the clean 
limit. 



B. Relaxation times 

The analysis of relaxation times is numerically diffi- 
cult. It turns out that the dominant relaxation times are 
very large near the first order phase transition. Accord- 
ingly thermal averaging is slow and measurements of the 
time decay of the density-density correlator have large 
statistical errors. However, although our data on relax- 
ation times is of rather poor quality, it is still of interest 
for comparison with the analysis of the vortex dynamics 
in the 2D system as well as for comparison with non- 
equilibrium measurements in YBCO. 

As in our previous analysis of the 2D system (35|, we 
measure the relaxation time scales in the layered sys- 
tem from the decay of the density-density correlator Cd 
from Eq. |To| in its Fourier transformed time-dependent 
form, Cd{k a b,k c ,t). We observe for high temperatures 
linear exponential decay of this correlator to zero for all 
k, where the k a b dependence of the decay time scales re- 
flects, like in 2D the hexagonal order in the system. 
The time scales decrease monotonically with k c for all 
k a b. For low temperatures, the time scales over which 
we measure the decay are often much smaller than the 
longest decay time scales themselves, so that decay of 
the correlator is only observable over a fraction of its ini- 
tial value. However, we can still, knowing that the vortex 
matter is liquid, extract time scales by fitting the data as- 
suming linear exponential decay. The longest time scales 
in the system are in all cases given by the decay of the 
3D Fourier component of Cd at the first reciprocal lattice 
vector in the a6-plane and fc=0 along the c-axis. 

Figure |ll| shows a plot of the dominant time scales de- 
pending on |q!t| 3 ^ 2 for different |a2T^|- The plots suggest 
very fast growth behavior, roughly 

r ~ exp(ci exp(c 2 |a T | 3/ ' 2 )), 

where c\ and ci are appropriate constants. We expect 
an exponential growth of the dominant time scales with 
the range of correlations, r oc expFl^ b , if the dynam- 
ics is activated, as it is in the 2D limit |39|. The data 
in Fig. ^ together with the low temperature relation 
ln(7 c ) tx 21nZ a ;, cx larl 3 ^ 2 [M, which is approximately 



also valid for our numerical data (see Fig. ||), strongly 
suggests such activated dynamics in the layered system. 

Figure |ll| shows evidence that time scales increase dis- 
continuously across the transition. If activated scaling of 
the type r ~ exp(ci exp(c2|ctr| 3 ^ 2 )) holds, the disconti- 
nuity in length scales should be amplified in the time scale 
discontinuity. The increase in t across the transition for 
|a2T?7|=0.05 is roughly a factor 500. In the case |a2T»y|=l 
the time scales increase across the transition according to 
our fits only by a factor 4. Although we expect the dis- 
continuity in time scales to decrease as |cv2T?7| increases 
and the end-point is approached, it is very likely that 
we strongly underestimate the time scales of exponential 
decay below the transition for |a2T??|=l, because we fit 
to a very small regime of decay (see inset of Fig. [ll]). As 
in the 2D case the relaxation behavior at very early 
times is faster than the final, linear exponential decay. 

As in the case of length scales, the data for all |a2T^| 
should collapse at low ax, i.e. in the continuum limit. 
This implies very fast growth of time scales below the 
transition ay for continuous systems with |a2T??| beyond 
the end-point. This extremely fast growth is reflected 
in the behavior of the system upon heating and cooling 
for |a2T7?|=3 (see Fig. 0), where strong hysteresis due 
to very slow decay of thermal fluctuations is observable 
below ax = — 8. 
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FIG. 11. Logarithm of logarithm of relaxation times plot- 
ted against «t| 3,/2 . Dotted and dot-dashed lines mark the 
width of the first order phase transition from the magnetiza- 
tion discontinuity for |q2t?7|=0.05 and |ck2t7/[=1 respectively. 
The inset shows logarithmic plots of original relaxation data 
with linear fits. 



1. Comparison with resistivity data 

Although the relevant time scales for dc resistance 
measurements in YBCO are non-equilibrium pinning 
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time scales which are not directly comparable to the re- 
laxation times measured in the simulation, some connec- 
tions can be made. While our time scales are irrelevant 
for flux flow conductivity, they are relevant in a highly 
viscous vortex liquid regime, where the time scales of 
thermal fluctuations or plastic deformations of the vor- 
tex structure are much larger then the pinning time. In 
this case one speaks about a pinned vortex liquid S . In 
such a liquid flux flow is strongly suppressed by pinning. 
When there is a large discontinuity of relaxation time 
scales at the first order transition, it is conceivable that 
the flux liquid becomes suddenly pinned and discontinu- 
ities in the resistive behavior occur. 

A decrease of the discontinuity in time scales as one ap- 
proaches the critical end-point agrees qualitatively well 
with the decreasing discontinuity in the resistive data as 
Bi c is approached [^6[ . The decrease and complete dis- 
appearance of the resistivity jump are difficult to explain 
with the assumption that the first order vortex lattice 
melting changes to a vortex glass melting transition at 
Bi c . A change in resistive behavior and a discontinuity 
in the resistivity at low currents should remain for such 
a transition || . 

In experiment a very steep drop of resistivity to zero 
with decreasing temperature is observed even for B < 
Bic, when there is no first order transition. This feature 
can in our picture be attributed to the doubly exponen- 
tially fast divergence of decay time scales with decreasing 
temperature once the extrapolation of the transition line 
is crossed. Although there is no time scale discontinuity, 
the system does a fast crossover from flux flow conduc- 
tivity to a strongly pinned liquid regime. 



2. Irreversibility and magnetization measurements 

A puzzling feature of the magnetization measurements 
near the first order transition in YBCO is the occur- 
rence of reversible-looking magnetization jumps without 
the corresponding latent heat |Tl[| which must accord- 
ing to the Clausius-Clapeyron relation (Eq. |J) occur for 
an equilibrium measurement. In Ref. pi large magneti- 
zation jumps which would correspond to entropy jumps 
of up to 11 Ais/vortex/layer, more than 10 times higher 
than the largest jumps observed in calorimetric measure- 
ments of approximately 0.8 fcs /vortex/layer are 
found without any evidence of a latent heat. 

If magnetization and specific heat data do not obey 
the Clausius-Clapeyron relation, this must be due to lack 
of equilibration in the system. The connection between 
seemingly reversible features in the magnetization data 
and irreversibility has in some cases been shown experi- 
mentally. For the vortex transition in BSCCO Farrell et 
al. |l4j found that the magnitude of the magnetization 
jump is for standard SQUID measurement techniques 
strongly correlated to irreversibility, even if irreversibility 



only becomes obvious at temperatures below the magne- 
tization jump. A similar effect in YBCO has been re- 
ported by Schilling et al. jl5| , where a change in slope of 
the seemingly reversible magnetization is shown to be an 
effect of irreversibility. 

Very recent data in YBCO by Ishida et al. in the sec- 
ond of Ref. p6[ shows magnetization jumps in a geometry 
B || a at B — 1.5T, while in this geometry the specific heat 
jump has in similar samples been observed to disappear 
between 6T and 4T [[ll| , a value in good agreement with 
our simulation data, if anisotropic scaling is assumed to 
apply. In Ref. Q| the ac (/=390 Hz) susceptibility is 
also measured and found to be an almost perfect image 
of the dc magnetization jump. This analogous behav- 
ior of the jumps in dc and ac susceptibility, of which the 
latter is clearly not a thermodynamic but a dynamic phe- 
nomenon due to non-equilibrium pinning effects fhSf , is 
suggestive of irreversibility effects even in the dc case. 

Although our simulation does not reproduce these ef- 
fects, the relaxation time scales we measure give a clear 
indication why such irreversibility related effects can oc- 
cur just below the extrapolation of the first order tran- 
sition line, i.e. at ut < —8. Let us assume that the 
growth in relaxation time scales is of the type r ~ 
exp(ci exp(c2|ar| 3 ^ 2 )) and further take an estimate of the 
slope of lnlnr below cvt = — 8 from the low temperature 



data in Fig. 11 for coupling |a2T??|=3.2 (beyond the end- 
point). Using the estimate 9(lnlnr)/9(|a T | 3 ^ 2 ) ~ 0.09 
let us consider a decrease in temperature from ax = — 8 
to ax = — 10. This corresponds to a temperature decrease 
of less than 0.3K below the the extrapolation of the first 
order line at a field of B \\ c = 0.25T or, according to 
anisotropic scaling with j — 7, to a field B || a of 1.75T. 
The considered temperature interval is thus comparable 
to the width of the magnetization discontinuity in the 
data from Ref. Gq|. For this decrease in ax we can ex- 
trapolate an increase in time scales from lnlnr sa 2.3 to 
lnlnr as 3.1, which is by a factor 10 5 . For an increase of 
time scales by a further factor 10 5 a further decrease of 
temperature by only 0.1 A" is needed. This implies that 
cooling by only fractions of a Kelvin below ay = —8 
for fields B < Bi c can cause the system to fall out of 
equilibrium. 

If one believes that falling out of equilibrium can 
lead to spurious magnetization jumps, then one may ask 
why a magnetization discontinuity is not always present, 
but in many cases disappears at low fields. For the 
data of Schilling et al. 
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|13J ] magnetization and specific 
heat data are found to be in good agreement with the 
Clausius-Clapeyron relation. A possible answer may be 
that the size of the magnetization discontinuity caused by 
non-equilibrium effects depends sensitively on the pres- 
ence of sample disorder. Such effects may make it negli- 
gible against the thermodynamic magnetization discon- 
tinuity at the first order transition in some samples like 
the ones used by Schilling et al. fl^,0, but not in others. 
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V. EFFECTS OF DISORDER 

We have investigated the effects of quenched random 
disorder on the layered system. We have limited our 
study to point disorder, i.e. the disorder realizations were 
always uncorrelated in different layers. We shall describe 
the effects of disorder on the first order transition as well 
as on correlations. 



A. Effects on the first order transition 
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FIG. 12. Plots of the order parameter density p upon heat- 
ing and cooling for different values of disorder. Solid lines 
(identical for all D) are for easier comparison. The hys- 
teresis decreases with increasing disorder, slightly faster for 
high |a2T?7|. For |a2T??| = 0.14 N c = 20. ..16 (increasing 
with D), for \a 2T v\ = 0.8 7V C = 60. ..40 (decreasing with D), 
and for |o2T»7| = 1.5 N c = 60. ..50 (decreasing with D). For 
\a 2T v\ = 1-5 p is offset by -0.025. 

Figure ^ shows plots of the order parameter density p 
upon heating and cooling for different values of disorder 
strength and coupling values varying from low coupling 
(|a2T?y| = 0.14) to values near the end-point (|a2T??| = 
1.5). We observe that the discontinuity persists in weak 
disorder, but decreases in amplitude and increases in 
width upon increasing disorder until it disappears en- 
tirely. The statistical noise makes it difficult to deter- 
mine the exact value of disorder for which the disconti- 
nuity disappears. However, it looks as if for D = 0.18 
the discontinuity near the end-point flo^T^I = 1.5) dis- 
appears, while for lower coupling values it is still present. 
This behavior would suggest that in a weakly disordered 
sample Bi c can be higher than in the clean case; for the 
case D = 0.18 we would expect B\ c between |a2T??| = 1.5 



and 0.8, i.e. roughly between 2 and 5 Tesla in YBCO. 

Our simulation does not capture the existence of an up- 
per critical field B uc . This can be interpreted as a sign 
that the appearance of B uc is strongly linked to the pres- 
ence of disorder clusters of a typical size, which are not 
present in a simulation with random Gaussian disorder. 
There is some experimental evidence for the link of B uc 
to the clustering of oxygen vacancies. On the one hand 
it has been shown that the location of B uc in YBCO is 
strongly correlated with the field strength for which the 
fishtail effect occurs at lower temperatures |l^]. On the 
other hand, the existence of the fishtail effect itself has 
been shown to be due to the presence of clusters of oxy- 
gen vacancies. A recent experiment by Erb et al. p7|l 
demonstrates that for constant oxygen vacancy density 
the fishtail effect can be suppressed by the destruction of 
vacancy clusters. 

This picture can be supported by the following qualita- 
tive argument. Let us assume the existence of a disorder 
length scale Id and consider the relation of Id to l m , the 
length scale which characterises the average vortex spac- 
ing and in the LLL also the order parameter coherence 
length. If l m > Id, the system is only sensitive to a spa- 
tial disorder average, and the disorder can be considered 
as effectively weak. When l m < Id, the system is sensi- 
tive to the length scale of disorder, and disorder effects 
are not weakened. We can define a field B uc ~ &o/l% 
Then the disorder is effectively weak for B < B uc , which 
is equivalent to the relation l m > Id- In our simulation, 
disorder is equally strong on all length scales, including 
length scales Id > l m an d may have effects similar to 
disorder clusters with a linear extent Id- The disappear- 
ance of the first order behavior in our simulation may 
thus be qualitatively comparable to the disappearance of 
the first order behavior observed in experiments at fields 
B > B uc . 



B. Effects on order parameter correlations 



By investigation of the order parameter correlations 
in the layered system near and below values of olt for 
which a first order transition occurs in the clean system, 
we may gain some qualitative information about the low 
temperature vortex state where disorder has destroyed 
the transition, i.e. for Bi CyCiean < B < B tc , (Bi CyCiean is 
the value of £>/ c in the clean limit) or for B > B uc . We 
choose to investigate the disordered system for a low cou- 
pling value, |a2T??| = 0.05, for which we see very large 
discontinuities at the transition in the clean system (see 
Fig. |and[ll]). 
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FIG. 13. Plots of the order parameter correlations with dis- 
order for \a 2T v\ = 0.05, D = 0.72, N ab = 72, iV c = 30: (a) 
vortex glass correlator, (b) phase correlations along the c-axis, 
(c) structure factor, (d) density correlations along the c-axis, 
(e) slowest time decay. All solid line fits are for ar = —8.35, 
and the corresponding correlation lengths are indicated. In 
(a) and (c) the dashed lines are fits for qt = —7.06. 

To probe if the system is in a vortex glass state we mea- 
sure the vortex glass correlator Q defined as the disorder 
average over 



C vg (r') = 



(^(rMr + r 7 )) (V>(r)V»»(r + Q) 
Q 4 



(13) 



where Q is denned in Sec. 0. C vg is a measure of the 
phase correlations in the system. We can determine the 
Fourier transform of this correlator from averages over 
the LLL coefficients as described in the Appendix, Sec. 



B 1 . C vg is in a clean system trivially short-ranged. It 
is gauge dependent but has typically a Gaussian decay 



behavior, m(C„ s ) 



18 , which can for D = be 



observed in our simulation. In the presence of disorder, 
decay with a finite, temperature dependent length scale is 
expected for a disordered vortex liquid Exponential 
decay with a correlation length l vg in real space would 
appear as 



1 



+ l V g 



(14) 



In Fig. [13] we show measurements of different density 
and phase correlators at a disorder strength for which all 
signs of a first order phase transition have disappeared 
from the heating and cooling measurements (not shown) . 
In Fig. [l^(a) l/C vg is plotted against k 2 . The behavior 
is in the plotted regime linear for (kl m ) 2 > 0.25. The 
data for wave vectors (kl m ) 2 < 0.25 corresponds to cor- 
relations more than a third of the circumference of the 
sphere. It is strongly affected by finite size effects and is 
therefore not taken into account for linear fits according 
to Eq. [ll|. The linear fits in Fig. |l3|(a) have rather large 
errors. Firstly insufficient averaging over disorder config- 
urations causes statistical noise. Secondly and more im- 
portantly the intercept — 1~ 2 depends sensitively on the 
exact choice of the linear fitting regime. It is however 
visible that for any choice it is finite and decreases with 
increasing \cct\- Although l vg grows with decreasing tem- 
perature, it is in the measured temperature regime only 
of the order of a lattice spacing. No long-range phase 
correlations typical for a vortex glass 0] are visible. 

The correlators we plot in Fig. |l^(b), (c), (d) and (e) 
have been introduced in Sec. IV. In Fig. |l3| filled sym- 
bols correspond to temperatures just above and at the 
location of the first order transition in the clean system 
and empty symbols to temperatures below it. In none 
of the static correlators is there any rapid change of be- 
havior across the location of the first order transition in 
the clean limit. However, in the relaxation behavior in 
Fig. [l3](e) an accelerated slowing down reminiscent of the 
discontinuity in the clean limit is still visible. 

All correlation lengths increase with \ax\- The corre- 
lation lengths l a b — 2<5 -1 Z m and l c from density as well as 
from phase correlations are reduced, for low temperatures 
by up to an order of magnitude, compared to the clean 
system (see Fig. ^|). Unlike in the clean system, Z c-p , i.e. 
the decay length of phase correlations along the c-axis, is 
distinctly larger than l Cj( i, the density-density correlation 
length. The phase correlation length in the plane from 
a fit of the vortex glass correlator at ar = —8.35 is of 
the same order as the density-density correlation length 
at the same temperature. Although olt = —8.35 is dis- 
tinctly below the phase boundary in the clean system, 
no signs of long-range phase correlations expected for a 
vortex glass phase are visible. 

We have not attempted a quantitative analysis of re- 
laxation behavior in the layered system plotted in Fig. 
|l3|(e), because we have not monitored the decay behav- 
ior over long enough time time intervals to see the domi- 
nant decay rates. However, we expect a behavior which is 
qualitatively similar to the the 2D limit of the disordered 
vortex liquid. We give in the following an analysis of this 
relaxation behavior, which extends and partly supersedes 
our earlier analysis given in Ref. pq |. We include with 
Fig. [l6| a quantitative correction to data from Ref. pjj] 
which was affected by an error in the simulation code. 
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FIG. 14. Typical relaxation behavior with increasing disor- 
der at Mm = 2.5 (near the first reciprocal lattice vector) and 
for an arbitrary different value of k. N a t = 72 and ol-it — — 9. 

Examples of the relaxation behavior of Cd(k,t) with 
and without disorder are shown in Fig. |l4| for the 2D 
limit. The relaxation behavior in samples with disor- 
der is richer than our earlier analysis |35j suggests. A 
seemingly faster decay than in the clean sample is visible 
at short times. However, this faster decay characterises 
only the initial decay behavior, which becomes clear in 
the data in Fig. |l4] for the case of strong disorder. At late 
times slow exponential decay on a time scale t/ due to 
pinning to the disorder becomes visible. To distinguish 
the two regimes of decay, we define an initial relaxation 
time Ti n by Ti n = Tin/ m(2), where Tin marks the time 
when the correlator has decayed to half its initial value. 
The final time scale tj can be found by exponential fits 
to the decay at late times. 
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FIG. 15. Initial and final relaxation time scales of the de- 
cay of the density-density correlator in reciprocal space with 
and without disorder. For D = system size N a t = 200, for 
D = 1.62 N ab = 72. a 2T = -9. 

Figure [l^ shows a typical logarithmic plot of the k de- 
pendence of initial and final relaxation times in the case 
of disorder compared to the clean case at the same tem- 
perature. The initial relaxation times Tj„ in the disor- 
dered case show qualitatively the same k dependence re- 



flecting the hexagonal structure as in clean samples [B5| , 
while the final relaxation times r/ do not depend on k. 
This seems natural, because the disorder contributions 
are on average equally strong for all k. 
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FIG. 16. (a) The inverse width of the first peak in the 
structure factor with linear fits for no and weak disorder in 
the low temperature regime. For olit = —10 and D — 
N a b = 200, otherwise N a b = 72. (b) Data of length and time 
scales plotted as ln(ri„) versus With linear fits for D — 
and D = 0.18. 

The dependence of the crystalline correlation length 
on Q.2T is plotted for different disorder strengths in Fig. 
[H^(a) and the dependence of the initial relaxation times 
on this length in Fig. |l6|(b). The initial relaxation times 
in the disordered system follow like in the clean system 
activated dynamics 35 1, in which the dominant relax- 
ation time scale grows linearly with the crystalline cor- 
relation length. 

The relaxation behavior in the disordered layered sys- 
tem plotted in Fig. looks similar to the 2D behavior 
at early times; the initial relaxation is faster than in the 
clean case. A more quantitative analysis of the relaxation 
behavior in layered systems would be of some interest in 
comparison to resistivity data in disordered samples. An 
analysis of the resistive behavior in YBCO before and 
after artificial introduction of point defects by electron 
irradiation jl9| concluded that transport in the disor- 
dered vortex liquid is dominated by viscosity rather than 
single vortex pinning, i.e. that the equilibrium relaxation 
time scales in simulations without a driving field could be 
relevant time scales. In the resistive measurements after 
irradiation in Ref. Jl9| no signs of a continuous vortex 
glass transition replacing the first order behavior before 
irradiation have been found. This absence of any thermo- 
dynamic phase transition in the case of strong disorder 
is in agreement with our simulation results. 



VI. DISCUSSION AND SUMMARY 

In our numerical work we have found a first order tran- 
sition which on the one hand coincides with the first or- 
der transition in YBCO, but on the other hand is not 
associated with vortex lattice melting. An analysis of 
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what triggers the decoupling transition at the particu- 
lar ax where it is observed is beyond the scope of this 
paper. The answer to this question could lie in the ob- 
servation made by Pierson and Vails |37]] that the first 
order transition in YBCO coincides with the onset of 3D 
LLL fluctuations. Because even strongly layered systems 
are 3D in nature below the first order transition, the on- 
set of 3D LLL fluctuations can be expected to occur at 
a fairly constant value of (Xt, and this may be reflected 
in approximately constant <xt along the phase transition 
line. 

Some readers may believe that the nature of the first 
order transition in YBCO is vortex lattice melting and 
dismiss the absence of a vortex lattice state below the 
transition in our simulation as an artifact of the spher- 
ical boundary conditions we use. In this case however 
it would seem astonishing that this first order transition 
associated with vortex lattice melting should quite hap- 
pily persist in our simulation in the absence of a vortex 
lattice. Besides, a first order vortex decoupling transi- 
tion [27] and another similar |30f thermodynamic transi- 
tion [ 24 1 , both not associated with vortex lattice melting, 
have also been observed in simulations of different models 
using periodic boundary conditions. 



A. Comparison with previous LLL-LD simulations 

From previous simulations using the same model, but 
with periodic boundary conditions, a single vortex lattice 
melting transition is reported p9pC| ] . We find there is a 
disagreement in the location of the transition between our 
simulation and Ref. |3Cj for weak couplings. This is cer- 
tainly due to the different choice of boundary conditions, 
because in the the 2D case the two choices of boun dary 
conditions are known to yield different results J32|,|36|,|31| . 
As the coupling increases, the disagreement in the loca- 
tion of the transition between our results and those from 
Ref. Q disappears. For couplings high enough to see 
the critical end-point, the LLL-LD model has never been 
investigated using periodic boundary conditions. 

The LLL-LD model simulations using periodic bound- 
ary conditions exhibit a vortex lattice state below the 
transition. If the transition was of the same nature for pe- 
riodic boundary conditions as for spherical layers, the ap- 
parent vortex lattice state could be accounted for by the 
use of system sizes smaller than the correlation lengths. 
(The largest system sizes in these simulations were of the 
order of only 40 vortices x 20 layers.) Small system sizes 
together with the restrictions in degrees of freedom im- 
posed by periodic boundary conditions may make a vor- 
tex liquid with very long length scales indistinguishable 
from a vortex lattice. 



B. Relevance for BSCCO 

We have so far only compared our numerical data 
with experiments in YBCO. Near the phase transition 
in BSCCO neither the LLL approximation nor the LD 
model are valid. The phase diagram of BSCCO has a first 
order transition line [p2|, which occurs at much smaller 
applied magnetic fields than in YBCO. Near the BSCCO 
transition H/H C 2 is of the order 10~ 3 . This means that 
the LLL approximation we use in our simulation cannot 
be expected to apply. The strongly anisotropic charac- 
ter of BSCCO is such that Josephson coupling between 
the layers may be negligible compared to electromagnetic 
coupling effects Q and so BSCCO is not accurately de- 
scribed by the LD model. 

However, some qualitative points of comparison can be 
made. It is for example noteworthy that in BSCCO the 
material parameters k as well as the layer periodicity d, 
T c and dB C 2/dT are of the same order as in YBCO, but 
typical estimates of the mass anisotropy 7 are between 
one and two orders of magnitude larger. This means that 
the end-point of our numerical transition line translates 
via the relation Bi c cx I/7 4 to fields that are between 
4 and 8 orders of magnitude lower than in YBCO. An 
experimental observation of a lower critical end-point is 
therefore not to be expected in BSCCO. 

In the low temperature limit, where all length scales 
are large, both YBCO and BSCCO should show the same 
universal behavior. Should our phase diagram be valid 
so that there is no thermodynamic vortex lattice melt- 
ing transition at a finite temperature below the experi- 
mentally observed first order transition in YBCO, then 
BSCCO should also be in a vortex liquid state below the 
first order transition. The consequence that the first or- 
der transition in BSCCO is not of a genuine vortex lattice 
melting character is in agreement with recent experimen- 
tal evidence that hexagonal neutron diffraction patterns, 
which signify that a vortex lattice or a vortex liquid with 
very long length scales, can be observed above as well as 
below the first order transition line fl54fl . 



C. Summary 

We have numerically calculated the phase diagram of a 
layered superconductor and found a first order transition 
line between two vortex liquid states. The length scales 
of order parameter correlations parallel and perpendicu- 
lar to the magnetic field as well as the longest relaxation 
time scales increase discontinuously at the transition but 
remain finite as the temperature is lowered. As the cou- 
pling strength between layers increases, the discontinu- 
ities in length and time scales decrease until the transi- 
tion line ends at a critical end-point. Shape, location and 
anisotropic scaling properties of the transition line and 
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its end-point as well as the size of the reversible mag- 
netization discontinuity are in excellent agreement with 
the experimental first order transition line and its low 
field end-point in YBCO. The approximate location of 
the end-point can be predicted from a qualitative argu- 
ment assuming that the transition is of a layer decoupling 
nature. However, the exact quantitative agreement of the 
location of the end-point with experiments by Schilling 
et al. JhJ could be a chance result due to cancellation of 
inaccuracies of our numerical model at low fields and/or 
finite size effects with disorder effects in real samples. 

Our results suggest that the transition in YBCO, which 
is commonly interpreted as vortex lattice melting, is of a 
liquid-liquid nature, with a low temperature vortex liquid 
phase in which length scales grow exponentially fast and 
time scales due to activated dynamics doubly exponen- 
tially fast with decreasing temperature. We argue that 
because not far below the first order transition the vor- 
tex liquid is highly viscous and effectively crystalline over 
large length scales, our picture can account for many ex- 
perimental features including resistance drops and Bragg 
peaks, which have so far been taken as evidence for the 
vortex lattice melting scenario. 

We have investigated the effect of quenched random 
point disorder on the system and found that the first 
order transition can be suppressed by strong disorder 
and that the low field critical end-point can be shifted 
to higher fields by weak disorder. We find no evidence 
of a thermodynamic transition to a Bragg glass phase or 
vortex glass phase replacing the first order transition line 
in the presence of strong disorder. 
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APPENDIX A: DETAILS OF THE 
HAMILTONIAN 



1. Quartic energy term 



with U p = Er=max(o, P -AO f( m ^P~ m) v m v p - m , with 
f(m,n) = A m A n (B(m + n + 1, 2N - m - n + 1)) 1/2 , 
where A n as defined in Eq. [I] and B is the beta function. 



2. Coupling to disorder 

Let the Gaussian random disorder O(r) in each layer 
be expanded as 



7 rp °" ' 



2=0 m=-l 



where the Y™ are spherical harmonics normalized with 
respect to the spherical surface of radius sj N/2. For 



real O, af 1 * 



a t m . The af 1 are chosen independently 



from a random Gaussian distribution with zero mean and 
variance |a™| 2 = D. Calculating the disorder correla- 
tions 0(r)0(r') and inserting A as defined in Eq. [?] fixes 
D = (rf Q 2 /fcsT) 2 A. For the integration of (layer indices 
omitted) 

iS^^f / d2r9(rMr)|2 (A2) 
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we define for < m < I and < q < N - 



I™q = / d 2 rY r 



(A3) 



, , (2l+l)(l+m)\ (-l) m . 

= A ^ f 2 J ( l m) ! [ -^rB(N- q+ i >q+m+ i) 

x 3 F 2 (m-l, m+l+1, q+m + 1; m + 1, N + m + 2: 1) , 

where 3-F2 is a generalized hypergeometric function. Sub- 
stituting these integrals into Eq. A2, truncating the 



contributions of noise components with frequencies / > 
( max = N (for a standard system size of N = 72 these 
are more than 10 20 times smaller than the largest con- 
tributing terms), and defining 
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The quartic coupling term in Eq. |^ can in each layer 
be expressed in terms of the LLL eigenf unctions in the 
following way (we omit the layer numbering indices): 



H 



quartic 



k B T 



^d [d*r\Tpf = ±r\U„ 



2N^'~ pl ' 

p=0 



(Al) 



yields 



n 



N N 

^ = EEv; 

q=0p=q 



t'o 



X c.c. 



(A4) 



Note that once g QiP is calculated the original coefficients 
a™ can be discarded and nc 
neously be held in memory. 



a™ can be discarded and need therefore not all simulta- 
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APPENDIX B: CORRELATIONS 

To compute correlations in reciprocal space, we per- 
form the spherical equivalent of a Fourier transform, the 
expansion in the discrete set of normalized spherical har- 
monics, YJ m (r). To a value of I corresponds k = l/R. 
Because the liquid is isotropic, the correlator in k space 
depends only on the magnitude of k, i.e. only on I and not 
on m. For better averaging we calculate the correlator 
for all m and average over the different m. 



1. Vortex glass correlations 

We calculate the Fourier transform of the correlator as 
defined in Eq. [l^ in terms of thermal averages of the LLL 
coefficients. 
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where the 7,™ are defined in Eq. A3 
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I N c lmin(N ab ,N ab -m) 



m=— !p=l \ q=max(0, — m) 
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1>p+n,ij'+m(i ) V p+n,q'(t )h,q' 
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where c signifies the connected average and the IJ^ are 
defined in Eq. A3. The high temperature limit of this 
correlator is for t = n — easily calculated analytically 
M as 



a^ C <W* >Q= iN- + 1 + 1)1 



(B3) 



For analysis of relaxation times this correlator is nu- 
merically Fourier transformed in the c-direction. To 
make C d (l/R,nd,t) converge in the continuum limit we 
choose £|| as unit of length when integrating along the 
c-axis: 

1 Nc 

C d (l/R,q,t) = — \^ CM/R,nd,t)cos(qxnd) (B4) 
where q takes values q = 2irm/ (N c d) for m = O...N c /2. 



2. Density-density correlations along the c-axis 

The real space density-density correlator for At = 
and Ar || c, Ar — nd is given by 



C d {nd) 



(\^(r)\ 2 \^{r + ndc)\ 2 ) 



f. 



Taking the spatial average over all r in one layer involves 
the same integrals over LLL eigenfunctions as the calcu- 
lation of the quartic energy term. Together with spatial 
averaging over different layers this yields 



3. Phase correlations along the c-axis 

The phase correlations along the c-axis as defined in 
Eq. [l^ are easily expressed in terms of the LLL coeffi- 
cients using the orthonormality of the LLL functions and 
the identity olt = 2iraH / f3Q 2 '- 
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where 



U n = 

P,Q 



(Bl) 



min(q,N ab ) 

^p^m^p+n^q—m > 

m— max(0,q— N ab ) 



with /(m, q — m) as defined in Sec. A 1 . 

The 2D density-density correlator in k space from ref- 
erence p2]] is easily generalized to three dimensions as 



APPENDIX C: FINITE SIZE EFFECTS 

The limitations in the region of parameter space for 
which we can run our simulation as well as the limita- 
tions in accuracy of our measurements are mainly due 
to limited availability of processor time. The simula- 
tion time grows, as seen in Sec. IV B, roughly doubly 



exponentially fast with decreasing olt- The ratio of cpu 
time to simulation time, given essentially by the num- 
ber of floating point operations needed for one update 
of the state, depends linearly on the number of layers 
N c , but due to the quartic energy term for large systems 
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quadratically on the number of vortices N a b. Finite size 
effects become important and therefore large system sizes 
necessary when correlation lengths grow large, which is 
unfortunately just in the regime of phase space where 
relaxation times also grow large, namely at low olt- 



1. Effects of limited N„ 



Most of our data have been taken using N a b = 72. For 
high temperatures finite size effects are negligible, be- 
cause there is little in-plane order and the range of crys- 
talline correlations in any one layer is much smaller than 
the spherical dimensions. The topological disorder due to 
the 12 topological defects imposed by spherical geometry 
p2| becomes at high temperatures negligible against the 
strong thermal disorder in the system. With decreasing 
a.T and growing l a b, finite size effects due to finite N a b 
become important. For any one ax, the effects of the 
limited number of vortices are the more severe the lower 
the coupling strength |a2Tf?| and the 2D effective temper- 
ature an are, because outside the continuum limit in- 
plane order increases with decreasing coupling strength, 
indicated by a smaller correlation length l a b (see Fig. |^) 
and smaller Pa (see Fig. ||). As a consequence of too 
small Nab, we see a decrease in length scales of crystalline 
correlations as previously observed in 2D systems J||] as 
well as a decrease of correlations along the c-axis. We in- 
terpret agreement of measurements taken with N a f, = 72 
and N a f, = 144 as indication that with N a b — 72 we are 
already close to the thermodynamic limit. 



2. Effects of limited 7V C 



Finite size effects due to limited numbers of layers are 
for equilibrium measurements not a severe problem, be- 
cause an increase of N c increases the cpu time only lin- 
early, and improves the rate of thermal averaging by the 
same amount. Therefore the disadvantage of increasing 
N c is only the increase in equilibration time as long as 
enough memory is available. If N c /2, the distance over 
which correlations decay, is not distinctly larger than the 
correlation length along the layers (see Sec. IV), we find 
that correlations both parallel and perpendicular to the 
c-axis are artificially enhanced and slower than exponen- 
tial decay behavior is observed, as visible in Fig. [Tt] for 
N c = 30. For equilibrium measurements of length and 
time scales, we always have N c > 10 x l c /d, which re- 
quired on occasion using N c up to 300 layers with 144 
vortices per layer. 



3. Thermal averaging and initial relaxation 

Initial relaxation has proved a very difficult problem 
in our 3D simulations, where relaxation times are so long 
that the entire simulation time is often limited to only 
few times the longest relaxation time. In principle, the 
relaxation can be started from an arbitrary state, where 
obvious choices are a random state or a ground state. 
The relaxation from a ground state has the disadvantage 
that it is very slow for low ar- On the other hand it is 
fairly easy to judge how far the system has relaxed when 
started with the same state in every layer. Because our 
system is always in a liquid state, a slower than expo- 
nential decay of correlations in the c-direction (which is 
in a sufficiently large system always removable by fur- 
ther equilibration) is a reliable sign of insufficient equi- 
libration. A case of insufficient relaxation can be seen 
in Fig. [ll] for A c = 180. When using systems started 
from a ground state, they are always allowed to relax 
for longer than the longest relaxation time r, in most 
cases longer than 5r. Relaxation from a random state 
has got the advantage of being faster than from a ground 
state. However, we found measurements from an insuffi- 
ciently equilibrated random state often indistinguishable 
from equilibrium measurements at a higher temperature, 
and therefore avoided starting equilibrium measurements 
from a random state. 

The most efficient way to obtain a well equilibrated 
system is to cool down or heat up a configuration ob- 
tained for a similar temperature and coupling strength 
and identical system size. We used this method when- 
ever such configurations were available. The cooling or 
heating rates in these cases are similar to those used in 
our hysteresis measurements. 



Cd(G) 




125 r/d 



FIG. 17. Typical measurements of static correlations along 
the c-axis. Open symbols are density correlations, filled sym- 
bols are phase correlations. Note for squares that correla- 
tions are enhanced by insufficient system size and for triangles 
the incomplete relaxation of the system. Statistical noise can 
for large r suggest slower than exponential decay (circles) as 
well as faster than exponential decay (stars). For all systems 
N ab = 72. 
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4. Loss of first order behavior 



The change of the hysteresis at the first order transition 
for decreasing numbers of layers can be seen in Fig. |l§| . 
The correlation length along the c-axis (see Sec. [V) just 
above the transition is for Io^t^I = 0.4 approximately 
6 layers. For N c = 20, correlations along the c-axis are 
enhanced significantly just above the transition. The sys- 
tem does not remain in the high entropy, decoupled state, 
and the hysteresis is lost. This effect is visible in Fig. [l8| 
in the reduction of the layer independence parameter T 
(defined in Sec. IIIB) just above the transition in case 
of the smaller system size. The loss of hysteresis with 
decreasing system size is rather sudden. Further increase 
in system size beyond 5-6 times the range of c-axis cor- 
relations affects the hysteresis measurements very little. 

In all cases the system sizes used for the hysteresis mea- 
surements are more than five times the correlation length 
in c-direction just above the transition. The correlation 
lengths are for selected values of |a2T?7| known from equi- 
librium measurements in large systems and otherwise es- 
timated by interpolation. We are aware that the system 
sizes used for the sweep measurements are in most cases 
not distinctly larger than the range of correlations below 
the transition. However, the observed size of the jumps 
does not seem much affected by this; the discontinuities 
in the order parameter magnitude seen in equilibrium 
measurements in much larger systems for selected |a2T??| 
agree well with the results from sweep measurements. 




FIG. 18. Order parameter density p and the measure of 
layer independence F upon heating and cooling. For insuffi- 
cient number of layers, F is reduced on the high temperature 
side of the transition and the first order behavior is lost. For 
N c — 20 p and T are offset by -0.05 (equal levels marked by 
solid lines). N ab = 72 and \a 2 TV\ = 0.4. 

A very important point to verify is that the loss of 
hysteresis at the critical point is not an effect of insuf- 
ficient numbers of layers. We have made sure that the 
system sizes used near the critical point are not only more 
than five times the correlation length, but also that the 
loss of the transition at |a2T??|=2.5 occurs in a system 
that is not only in terms of layers, but also in terms 
of the natural length scale larger than the system for 



|o;2T?7|=2, where a transition is still visible. For \a2TV\ = 2 
and N c d/£\\ = 65 we see a clear transition, while for 
|a2T??|=2.5 and N c d/£\\ = 80 there is no sign of a transi- 
tion in the range —8.3 < ax < —7.3. 
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FIG. 19. Order parameter density p and Abrikosov num- 
ber /3a upon heating and cooling. For insufficient number 
of vortices per layer, /3 a is increased on the low temperature 
side of the transition and the first order behavior is lost. For 
Nab = 66 and N a b = 54, p is offset by -0.05 and -0.1 respec- 
tively and Pa is offset by -0.015 and -0.03 respectively (equal 
levels marked by solid lines). N c — 12 and |«2T??| = 0.05. 

If the number of vortices per layer N a b is reduced, the 
hysteresis at the transition first decreases and then dis- 



appears. For the general reasons outlined in Sec. C 1 
the effect becomes stronger as |a2T?7| decreases. Below 
some critical N a b, which increases with decreasing I^t^Ij 
the in-plane order below the transition, reflected in the 
Abrikosov Ratio Pa, is so much affected by the spher- 
ical topology that the discontinuity in (3a vanishes and 
the transition disappears. This behavior is shown in Fig. 
|Tg| for |ct2T^| = 0.05. The size of the discontinuities de- 
creases noticeably between N a b = 72 and N a b = 66. For 
N a b = 54 the transition has disappeared. For N a b = 72, 
we still measure a transition at |aaT^| = 0.02, but not 
at |cf2T?7| = 0.01. Thus the need of increasing N a b limits 
the exploration of the phase diagram for very low |a2T??|- 

We estimate that the limitation of the number of vor- 
tices we study to N a b — 72 affects our hysteresis mea- 
surements for |c«2T?7| < 0.05. We can detect size depen- 
dence in the location of the phase transition for system 
sizes N ab = 72 and N ab = 144 only for \a 2 Tr]\ < 0.05 
(see the phase diagram in Fig. |^). The hysteresis mea- 
surements plotted in the same figure show that the size 
of the discontinuities at the phase transition does not 
change noticeably between N a b — 72 and N a b — 144 for 
|c*2T?7| = 0-14. For |«2T??| ^S> 0.14, which applies to the 
region near the critical point, the system should at the 
phase transition be well simulated using N a b — 72. This 
is confirmed by the agreement in the magnetization dis- 
continuity between N a b — 72 and N a b = 144 in Fig. ||(b) 
for T=89.3. 
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